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On the rise of small air bubbles in water 


By P. G. SAFFMAN 


Trinity College, Cambridge 


(Received 9 February 1956) 


SUMMARY 

his paper is concerned with the motion in water of air bubbles 
whose equivalent spherical radi are in the range 0-5—4-0 mm. 
These bubbles are not spherical but are, approximately, oblate 
spheroids; and they may rise steadily in a vertical straight line, or 
along a zig-zag path, or ina uniform spiral. ‘lhe rectilinear motion 
occurs when the radius is less than about 0-7 mm, and the other 
motions occur for larger bubbles. ‘There is disagreement in the 
literature as to whether it is the zig-zag or the spiral motion that 
occurs. It was found experimentally that, when the bubbles are 
produced in the manner described in this paper, only the zig-zag 
motion occurs when the radius of the bubble is less than about 
Imm, but bubbles of larger radius either zig-zag or spiral 
depending upon various factors. 

The spiralling bubble is treated theoretically by assuming that 
the flow near the front of the bubble is inviscid (the Reynolds 
number of the motion is several hundred) and considering the 
distribution of pressure over the front surface. Equations are 
obtained relating the geometrical parameters of the spiral, the 
shape of the bubble and the velocity of rise. ‘The analysis is 
simplified by assuming that the pitch of the spiral is large compared 
with its radius, and the velocity of rise and shape of the bubble are 
determined as functions of the radius. The experimental and 
theoretical values are compared, and fair agreement found. 
Reasons to account for the disagreement are proposed. 

A modification of the theory is proposed to take account of 
the presence of impurities or surface-active substances in the 
water, and the velocities of rise thus predicted are in agreement 
with the experimental observations. 

The zig-zag motion is treated in a similar way, and the analysis 
leads to an equation which determines the stability of the rectilinear 
motion. The value of the Weber number at which the rectilinear 
motion becomes unstable is deduced, and is found to be in fair 
agreement with experiment. ‘The experimental evidence on the 
wake behind solid bodies is described briefly, and reasons are 
given for suggesting that the zig-zag motion is due to an interaction 
between the instability of the rectilinear motion and a periodic 
oscillation of the wake. 


F.M. 
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1. INTRODUCTION 


Many authors have studied experimentally the motion of air bubbles 
in water and other liquids. Datta, Napier & Newitt (1950) and Peebles & 
Garber (1953) describe briefly many of these experiments, and further 
experiments have been carried out by Haberman & Morton (1953). It has 
been observed that when the bubble radius is about 0-7mm and the 
Reynolds number is several hundred, the character of the motion changes. 
Bubbles of radius less than 0-7 mm rise steadily in a straight line; but, 
when the radius exceeds 0-7 mm, the bubble either zig-zags from side to 
side or rises in a spiral which appears to be quite uniform. (It is convenient 
to specify the size of an air bubble by its equivalent spherical radius, this 
being the radius of a sphere with the same volume as the bubble. Whenever 
the radius of a bubble is mentioned in this paper, the equivalent spherical 
radius is implied.) ‘The various experimenters are not in agreement 
as to whether the bubbles spiral or zig-zag, and an experiment was 
performed (described below in $2) in order to investigate this matter 
further. 

The reported experimental evidence indicates that the frequency of the 
zig-zag motion is about seven per second and independent ot the size of the 
bubble; the angular velocity of the spiralling motion also appears to be 
independent of the size and is about 30 radians per second, which corre- 
sponds to about five revolutions per second about the axis of the spiral 
The radius of the spiral appears to decrease slowly as the bubble size 
increases, being about 1-5 mm when the radius of the bubble is about 
1:0 mm; the amplitude of the zig-zag is also of this order of magnitude. 
As the radius of the bubble increases past 2-0 mm, the spiral or zig-zag 
begins slowly to disappear; and by the time the bubble radius is about 
3-0 mm, the spiral or zig-zag motion has disappeared completely. The 
motion is then roughly rectilinear, and continues to be so as the 
bubble radius increases, although the shape of the bubble alters 
considerably. 

Small bubbles of radius less than 0-5 mm appear to be spherical. 
Larger bubbles are distorted, and the distortion becomes greater as the 
size increases further. Bubbles of radius between 0-5 mm and 3-0 mm 
are oblate spheroids approximately, but the shape of the larger bubbles 
in this size range tends to fluctuate somewhat. For a bubble of radius 
larger than 3-0 mm, the shape is very irregular and fluctuates considerably ; 
then, at a radius of about one centimetre, the bubble has the shape 
of a spherical cap, and this shape is maintained as the bubble radius 


increases. 

The velocities of rise have been measured for spiralling and zig-zaggin 
bubbles. ‘The measurements show a fair amount of scatter, but there 
appears to be no systematic difference between the velocities of spiralling 
or zig-zagging bubbles, i.e. the velocity appears to be independent of the 
nature of the motion. Several authors (e.g. Haberman & Morton 1953) 
have observed that the velocities of bubbles in the size range 0-3-3-0 mm 
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are reduced somewhat if the water is impure or contains ‘ surface-active ’ 
substances. However, Haberman & Morton report that air bubbles 
in filtered tap water rise with the same velocities as bubbles in distilled 
water. 

The theoretical work so far published has been confined to bubbles 
which are either so small or rise so slowly that the approximations of the 
slow viscous motion theory can be applied, or to large bubbles which are 
in the form of spherical caps (see Davies & Taylor 1950). It is the object 
of this paper to give a theoretical account of some aspects of the motion of 
the spiralling and zig-zagging bubbles. 

In $3 we consider the spiralling bubbles. ‘The method used is briefly 
as follows. We suppose that the motion of the water near the front of the 
bubble is inviscid and irrotational, and calculate it in terms of the velocity 
of rise, the shape and orientation of the bubble, and the dimensions of the 
spiral. We can then find the pressure at the front surface of the bubble. 
The pressure inside the bubble is the sum of the pressure outside and a 
term involving the surface tension and curvature of the surface. We may 
take this sum to be constant at all points of the surface, since the pressure 
variations inside the bubble can be neglected owing to the small density 
of air compared with that of water. We apply this boundary condition 
and obtain relations between the quantities mentioned above; in particular 
we obtain an equation giving the velocity of rise of a spiralling bubble in 
terms of its size. It is to be noted that this method does not require a 
knowledge of the boundary condition to be satisfied by the velocity of the 
fluid at the surface of the bubble, since it is assumed that the effect of 
viscosity on the motion at the front of the bubble is confined to a thin 
region. 

This method is used in §4 to investigate the motion of a zig-zagging 
bubble. In this case, we find that instead of obtaining an equation involving 
the velocity of rise and the bubble size, a stability equation results from the 
analysis. ‘This equation shows that the motion near the front of a bubble 
which is rising steadily in a straight line is unstable when the bubble is 
sufficiently oblate. 

To obtain further information about the motion, it is necessary to 
consider properties of the wake behind the bubble. ‘There 1s some relevant 
evidence on the wake behind solid bodies, and this is discussed in §4 in 
connection with the zig-zag motion. 


2. EXPERIMENTAL METHOD AND RESULTS 

Air bubbles were observed as they rose up through a rectangular tank 
of height 15 in. and cross-section 12 x 12 in., which was filled with either 
filtered (standard No. 1 filter papers were used) or unfiltered tap water at 
temperatures between 16° C and 19°C. They were released at the bottom 
of the tank from the end of a capillary tube, and were introduced into the 
other end of the tube by squeezing a rubber syringe containing air (see 
figure 1). 
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The experimental procedure was as follows. ‘lap ‘A’ was opened, 
and tap ‘ B’ was so adjusted that water flowed very slowly from the ‘l’-shaped 
glass tube through the capillary tube into the tank. ‘Then tap ‘A’ was 
closed, and the rubber syringe containing air was squeezed gently until an 
air bubble of about the size required formed in the funnel-shaped end of 
the capillary tube. Tap ‘A’ was opened again, and the air bubble 
rose up along the capillary tube until it came opposite to the scale. ‘Tap ‘ A’ 
was then closed; the bubble stopped moving, and its length was measured 
on the scale. Since the bore of the tube had been found accurately by the 
usual methods, the volume of the bubble could be calculated with sufficient 
accuracy. ‘lap ‘A’ was now opened again and kept open until the bubble 
reached the end of the capillary tube at the bottom of the tank, whereupon 
the tap was closed so that the water in the tank could come to rest. ‘Tap ‘ A’ 
was opened once again, and the bubble rose up through the tank and was 
observed either through the sides of the tank or from above. Once sufficient 
skill in making the delicate adjustments had been acquired, the method 
worked very well, and it was possible to obtain bubbles of about the desired 
size at will. With this arrangement it was possible, by adjusting the tap ‘ B’, 
to vary the pressure in the ‘I’-shaped tube in such a way that bubbles in the 
capillary tube could be moved, with adjustable speed, in either direction. 


FROM MAINS waTeQ 
od Tap 
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Figure 1. Experimental apparatus (not to scale). 


Bubbles of radii between 0-5 and 2:3 mm were observed during the 
course of the experiments, a capillary tube of bore 0-8 mm being used 
for the smaller bubbles and one of bore 1-6 mm for the larger. Bubbles 
of radius larger than 2-3 mm could not be produced since they always 
broke up before leaving the capillary tube, even when tubes of larger bore 
were used. In fact, great care was needed in order to make bubbles of 
radius larger than about 2:0 mm. The following observations were made 
on bubbles in filtered water. 

If the radius was less than 0-7 mm, the bubble rose steadily in a straight 
line. If the radius was between 0-7 and 1:0 mm, the motion was a zig-zag, 
i.e. a side-to-side movement in a plane; the direction of the plane remained 
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constant as the bubble was rising, but there was no apparent tendency for 
the bubbles to prefer any particular plane. 

When the radius was between 1-0 and 2:3 mm, the situation was more 
complicated. It was found that bubbles in this size range would either 
spiral or zig-zag, depending upon various factors. For instance, it sometimes 
happened that a large bubble would break up into two parts just before 
it left the capillary tube. ‘The first part would rise through the tank and 
always spiral, whilst the other stayed in the capillary tube. The size of 
the part left behind could be measured by bringing it back to the scale; 
and, since the size of the original bubble had been measured, the size of 
the spiralling bubble was known. When this break-up occurred, the radius 
of the bubble that left the tube was always greater than 1-3 mm; it was 
not possible to produce this break-up at will and smaller bubbles could 
not be produced in this manner. 

However, if care was taken to see that external disturbances were 
reduced to a minimum and the bubble left the capillary tube in one piece, 
then the bubble would usually zig-zag. A bubble could usually be made 
to spiral by ‘hitting’ it, 1.e. by placing an obstacle above the end of the 
capillary tube and in the path of the bubble. Further, if two bubbles were 
released shortly after one another, the second would spiral if the first one 
did; but if the first one zig-zagged, then so usually did the second. It was 
noticed that a zig-zag would sometimes change into a spiral, but the opposite 
was never observed. Bubbles of radius less than 1-0 mm were always 
observed to zig-zag, even after being ‘hit’ or released in the wake of a 
spiralling bubble. Bubbles were also projected at speed from the capillary 
tube, but there appeared to be no correlation between the initial velocity 
of the bubble and the occurrence of zig-zag or spiral motion. 

The general impression received from these observations was that, 
when the radius was greater than 1-0 mm, the zig-zag motion was unstable 
to sufficiently large disturbances and the bubble would spiral, the tendency 
to spiral increasing as the size increased. 

Furthermore, bubbles of radius about 0-7 mm were sometimes observed 
to rise steadily, then zig-zag, then go straight again, zig-zag again, and so 
on. When bubbles of this size were released in the wake of a spiralling or 
zig-zagging bubble, they would zig-zag with the usual frequency, but 
‘hitting’ these bubbles would not produce a regular zig-zag. Slightly 
smaller bubbles of radius about 0-6 mm were also ‘hit’ and released in 
the wake of a larger bubble, but they always rose steadily in a straight line. 

Similar experiments were carried out using unfiltered tap water. The 
results were the same except that the zig-zag motion first occurred at a 
radius of 0-8 mm, and that it was not possible to produce the spiralling 
motion when the radius was less than 1-1 mm. 

To sum up, these experiments suggest strongly that when the bubble 
radius reaches the value at which the steady rise in a straight line becomes 
unstable, the zig-zag motion appears first; and that at a slightly larger 
radius, this zig-zag motion is itself unstable to sufficiently large disturbances, 








254 P. G. Saffman 


and a spiralling motion appears. ‘lhe observations that the zig-zag motion 
occurs first and sometimes changes into the spiralling motion, whereas the 
spiralling motion never changes into the zig-zag motion, indicate that the 
spiralling motion arises from a later instability than the zig-zag motion. 
However, some authors (e.g. Miyagi 1925) claim to have observed only 
spiralling bubbles, and it seems likely that under certain experimental 
conditions, zig-zagging does not occur for any bubble size. ‘These conditions 
will presumably be related to the way in which the bubbles are produced, 
and it is possible that, in such cases, the disturbances that are present are 
sufficiently strong to make the zig-zag motion always unstable. 


3. THE MOTION OF A SPIRALLING BUBBLI 
The formulation of the problem 

We shall now give a theoretical treatment of the motion of a spiralling 
bubble. ‘The problem is complicated, and we shall make several simplifying 
assumptions in order to render it tractable. One complication is the need 
to determine the shape of the bubble. It can be shown that a bubble is 
nearly spherical, when the Reynolds number is large, if 2¢R/W? (the inverse 
of the Froude number) and 2pRW?/T (the Weber number) are both much 
less than one. Here, 7 is the surface tension, R the bubble radius, W the 
velocity of rise, and p the density of the fluid. For spiralling (and zig-zagging) 
bubbles 2gR/W7* is small but 20RW?/T is of the order of unity ; consequently, 
we cannot suppose the bubble to be nearly spherical and make use of the 
simplifications that are then possible. ‘This conclusion is consistent with 
the experimental observation (see, for example, Rosenberg 1950) that 
spiralling bubbles are not spherical but are approximately oblate spheroids. 
It will be assumed, for the purposes of the analysis, that the bubble is an 
oblate spheroid. ‘This is not strictly true, but it seems to be a reasonable 
approximation. 

The spiral along which the bubble moves appears from the experiments 
to be roughly uniform and to be described at a constant rate. We shall 
therefore suppose that the motion is steady when referred to axes fixed in 
the bubble, and that the shape of the bubble is constant. We can, in fact, 
imagine the bubble to be at rest and the water to be streaming past it steadily. 
In order to simplify the analysis further, it will also be supposed that the 
angle of the spiral (which is defined as the angle between the upward 
vertical and the tangent to the spiral) is small, and we shall neglect, in the 
subsequent equations, terms that are of the same order of magnitude as this 
angle. Later, we shall estimate the numerical magnitude of such terms 
and see that this approximation is reasonable. It will also be supposed that 
the angle between the axis of symmetry and velocity of the bubble is small 
compared with the angle of the spiral; and it then seems permissible, for 
the sake of simplicity, to take the axis of symmetry coplanar with the velocity 
of the bubble and the vertical, i.e. to suppose that the axis of symmetry 
lies in the tangential plane of the cylinder that contains the spiral. This 
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assumption is consistent with the rough observations made by several authors 
(e.g. Miyagi 1925) that the axis of symmetry and velocity of the bubble are 
parallel, and can also be shown to be self-consistent with the results of the 
analysis given below. 

We have now set up the geometry of the problem, and we shall now 
consider how we may calculate the pressure p in the water. ‘The basic 
assumption, upon which the analysis in this paper is based, is that the 
pressure over the surface of the bubble near the front stagnation point is 
the same as if the whole flow were inviscid and irrotational. In effect, we 
assume here that the motion near the front of the bubble is inviscid and 
that the velocity and pressure near the front stagnation point are primarily 
determined by the shape of the neighbouring surface. Now the boundary 
condition to be satisfied by the pressure is (neglecting pressure variation 
inside the bubble) 

p+7T(1/R,+1/R,) = const. (1) 


over the whole surface of the bubble, where R, and R, are the principal 
radii of curvature of the surface. We then make use of the above assumption 
by supposing that the first and second derivatives of equation (1), differen- 
tiating along the surface, are zero at the front stagnation point when the 
pressure calculated on the assumption of ideal flow is substituted for p. 
From the equations so obtained, we shall derive equations relating the 
geometrical parameters of the spiral, the shape of the bubble and the 
velocity of rise. 

It is worth noting that Davies & ‘l'aylor, making assumptions about 
the pressure that are basically the same as those made here, calculated the 
velocity of rise of large spherical-cap bubbles, and the theoretical velocities 
agreed well with the velocities which they observed experimentally. ‘The 
Reynolds number for these large bubbles was several thousand. ‘The 
Reynolds numbers of the spiralling bubbles dealt with here are several 
hundred, which should still be sufficiently large to justify the assumption 
that flow over the front of the bubble is essentially inviscid. We shall, 
however, modify later the assumption of ideal flow in order to take account 
of possible separation of the flow around the bubble. In this way, we shall 
obtain results in agreement with the observations of bubbles rising in 
impure water. ‘his point will be discussed again in more detail. 

It should be noted that the theoretical treatment given here is incomplete 
in the sense that we do not investigate conditions at the rear of the bubble, 
owing to the extreme complexity of the motion there. ‘he consequences 
of this incompleteness are that we shall obtain fewer equations than 
unknowns, and that we shall obtain no information as to the causes of 
the spiralling motion. 


The calculation of the pressure 
We denote the vertical velocity of rise by W, the angular velocity with 
which the spiral is described by ©, and the radius of the spiral by d. 








256 P. G. Saffman 


We take axes (x, vy, 2) with origin fixed at the centre of the bubble; x is 
the axis of symmetry, and the plane x, vy contains the vertical. 

Further, we take oblate spheroidal coordinates yu, ¢, w defined by 
x= kul, _ R(1 pe”)! (14 {*P “cOsS@, 2. = R(1 —p?)! (1 + (?)} 2sinw (see, 
for example, Lamb (1932) § 107). We suppose that the surface of the bubble 
corresponds to € = ¢, so that the eccentricity of a meridian plane of the 
bubble is (1+ ¢?)-!, and the equivalent spherical radius is k{,(1 + ¢?)}¥3. 
We denote by y the angle between the axis of symmetry of the bubble and 
the upward vertical. 

We now reduce the bubble to rest by superposing a velocity — W and 
an angular velocity — 2 on the fluid, where W = (W cos y, Wsinx, 0) and 
Q = (Qcosy, Qsiny, 0). The velocity of the undisturbed fluid is now 

W-2x(d+r), where d = (0, 0, d), and r = (x, y, 3). Let us denote 
by ® the velocity potential of the disturbance due to the bubble, so that 
V = VO—-W-Qx(d+r) is the velocity of the fluid relative to the axes 
(x, y, 3). (The actual velocity of the fluid relative to axes fixed in space 
is, of course, just V@.) 

The surface of the bubble is now a streamline, so that V.n = 0 on 
where n is the outward normal to the surface. Hence, after some 


S$ = So» 


algebra, we obtain 
c® — ; 
(=) = pk(W cos y + Qdsin x) + 
e : < , , 


+k(1—p?)2(1 + 2)-126,(W sin y — Qd cos x)cos w + 
+ OR? n(1 — p?)2(1 + G) PF sin ysinw. (2) 
We are supposing the fluid to be incompressible, so that ® satisfies 
V*@ = 0. Suitable solutions of this equation which fit the boundary 
condition (2) and which vanish at infinity are 
p(1—cotc-12),  (1—p?)¥8(1 + 0%)-"4[¢ — (1+ 2?) cot C]cosw, 
w(1 — pw?) (1 + C7) 2[3 6 cot €-3.+ (1+ 2?) ]Jsinw. 
We then find, after some algebra, that on € = G, 
M = —kp(1+f)2(1—C cot f)(W cos x + Qdsin y)Z-1 — 
—RZYVE(1+ &)-'*(Wsin yx — Qd cos y)(1 — p?)!* cos w + 
+ QR? sin y(C? + 1)!*(X —Qu(1—p?)!*sinw, (3) 
where for brevity we write 
3€cot ?(-34+(1+@)"? 
6C? + 3)cot?f-66-—f(11+ 2)! 





’ 


elie ial’ 


VY = (242-14 2)cor tg, 
Z = (1+2){(1+ 2)cot ff}, 
and the suthx of ¢) has been dropped. 
Let V,, and V, denote the components of V on the surface of the bubble 
in the directions of increasing 4 and w respectively; V, is, of course, zero. 
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We use (3) and obtain, after some algebra, 
(wm? + C?)P2V = —Z-(1+ C?)(Wcos x + Qdsin x)(1 — p?)!2 + 
+ 2Y(1+ ¢?)'2(Wsin x — Qd cos yx) cos w + 
+ Qksin y (1 + C?)#2[ —-6+(1—2p?)(X—2)] sinw, (4) 
and 
V, = 2Y(Wsin yx — Qd cos x) sin w — 
— Ok cos x (1+ C?)82(1 — p?)2 + OkRsin y Xucosw. (5) 


At the front stagnation point, V,, and V,, are zero, and the coordinates of 


the stagnation point, (,, w,) say, are given by putting the expressions (4) 
and (5) equal to zero. 
The pressure is then given by* 
p/p t+ V2 +2¢H — $Q?m" = const. 
on a streamline, where g is the acceleration due to gravity, 
= [(d+2z)?+(xsin y—ycos x)*}!? 
denotes the perpendicular distance from the axis of the spiral and 
H = xcosy+ysin x 
is the vertical height above the centre of the bubble. Further, 
Lot (e+e c 
R, RR, (482 RCC +I! 2) 


by well-known results of solid geometry. 





The application of the pressure boundary condition 


Now the surface of the bubble is a streamline and so, in view of the 
initial assumption about the pressure, the seine condition is 


2 Soe, een, on I\\ =~ 0. 6 
(5. Cw’ Op?’ Cucw’ x3) ( a ” i “” 2 ( ) 


the derivatives being evaluated at the stagnation ae aie w,). 





From (6), ¢/du gives 
TCpn(40? + p? + 3) 
pk(C? + 7)? bi sd 
=a O*kd(1 + c*) 2u(1 — wy 2 sin w 4 
+ OPk? (ul? sin?y — (1 + C?)(sin’w + cos*y cos’w) + 
+ (2u?—1)C(1 + ?)!? sin ycos x(1—p?)"? cosw} = 0; (7) 


—gkicos x + gksin x (1+ @)'2u(1 — p?)-?? cos w — 





“This equation is the appropriate form of Bernoulli’s equation for steady motion 
relative to moving and rotating axes. One way of deducing it is as follows. The 
pressure is given by p/p~+}(V®)?~@®/0t--gH = constant, where O/dt denotes 
differentiation with respect to axes fixed in space. The velocity in space of the axes 
(x, vy, 2) is W- Q»d and the angular velocity of these axes is 92; and it follows 
that @®/dt —[W- Q»(d- r)].V®. Remembering that V® = V+W--Qx (dr 


we obtain the given equation on substitution. 
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and d/dw gives 

h(1 + C?)¥#(1 — ?)¥4{ 9 sin y sin w + Q2d cos w + 
+ Ok sin?y (1 + ¢?)"2(1 — w?)"? sin w cos w + Qk sin x cos x Gusinw} = 0, (8) 
the suffixes of », and w, having been dropped since all expressions are 
evaluated at the stagnation point. 

We now suppose that powers of yx and (1—,?)!? higher than the first 
can be neglected. It follows from equations (12) and (4’) obtained below 
that these quantities are of the same order of magnitude as the angle of the 
spiral, which is supposed to be small. Equations (4), (5), (7) and (8) thus 


become 

















Wz-(1 + @7)?2(1 — p?)'2 + 2Y cosw (Wy —Qd)—-ORXysinw = 0, (4’) 
— OrR(1 + C?)47(1 — p?)!? + 2Y sinw (Wy — Qd)+ QRXycosw = 0,  (5’) 
4TC . QeR (+2 (, — “) 2 (2d <4 
Q2RC 
_ (1 se ) xcosw = UV, (7’) 
(2?7RkC (2d 
and (1 + a sin w +4 ~ cosw = 0). (8’) 
We now eliminate « and w from these four equations and obtain 
; : =) ( kgZyx okXy 
2 = SN a OEE 
= (u X, 4 Qd ’ (9) 
and 
WwW [42 f OAd\12 
7 (1+¢ x’ r =| 
_ Ay _ ATE, [Oa Maen 
me Qd- pgk?( 1+2) c) go X a od) | p 7) 


where we have neglected terms of order (2?k¢/g and 0?k?/W?, which can 
be shown to be of the same order of magnitude as the angle of the spiral. 
Let us now consider the second order derivatives in (6), using again the 
assumptions that y and (1 —,?)!? are small. From 0?/du? we obtain 
we 7 4TL 
— ca hie i . (1 1) 


gk «1+ | pgk®(1+ 2)? 
from ¢?/cucw we obtain a combination of equations (4’), (5’) and (8’); 
and from ¢?/¢w? we obtain a combination of equations (11) and (7’). 

Thus, we have obtained three complicated equations, namely (9), (10) 
and (11), which give the five unknowns W, Q, d, ¢, y in terms of T/p, g and 
the equivalent spherical radius r (ry and & are related by r? = k°%(1+ @)). 
These equations have been derived from an analysis of the motion at the 
front of the bubble; and it is not surprising that we should obtain fewer 
equations than unknowns since we have considered only one region of the 
flow. Other equations would presumably result from a study of the motion 
in the wake and the forces acting on the bubble. 

However, we can simplify these equations and determine ¢ and the 
velocity of rise provided that we neglect kgZy/(QdW). Later, we shall 
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make an estimate of this quantity and see that its neglect is reasonable for 
bubbles whose radius is less than 1-5 mm, but not for larger ones. 
Equations (9) and (10) then become 








2¥(W- = \= = hgX, (9") 
x X 
and 
W  gkXx 4T¢ 
Ser S ay 
Combining (10’) and (11), we obtain 
Qd xX 
RZ (4) 
and from (12), (9’) and (11), we then obtain 
T 1+@ 
; ao —" (13) 


pek(1+¢7)? = 2Y¢(Z-X) 

Equation (13) determines the shape parameter ¢ as a complicated 
function of the bubble radius, and then the velocity of rise follows from 
equation (11). Equation (12) contains the remaining three unknowns, and 
gives the ratio of the angles made with the vertical by the velocity and 
axis of symmetry of the bubble. 


Comparison with experiment 

We shall now discuss the above theoretical results and compare them 
with the experimental observations of various authors. It will be seen 
that the agreement between the velocity of rise and bubble shape predicted 
by equations (11) and (13) and the experimental observations is reasonable, 
considering the nature of the approximations that have been made. We 
shall then examine some of the approximations in more detail and see how 
they might account for some of the disagreement between theory and 
experiment. 

Figure 2 shows how X, Y-! and Z vary with ¢. It will be noticed that Z 
is greater than X only if € is less than 1-52. Hence, the right-hand side of 
equation (13) is positive, which is obviously necessary if 7 is to be real, 
only if the bubble is sufficiently flattened (it will be remembered that the 
smaller ¢, the less spherical is the bubble). In other words, no bubble can 
spiral unless it is flatter than the oblate spheroid corresponding to ¢ = 1-52. 

In figure 3 is shown the solution of equation (13) for air bubbles in water, 
where the typical values T/p = 74cm*sec~? and g = 981 cmsec~ have been 
used. For each value of r less than 3-6 mm, there are two values of ¢; 
and we take the larger value as the one corresponding to physical reality, 
since then ¢ decreases as r increases. (Surface tension tends to keep the 
bubble spherical, and the magnitude of the surface tension effect is 
inversely proportional to the radius of the bubble.) 

Further, there are no values of ¢ which correspond to r greater than 
3-6 mm, and we can interpret this as meaning that the spiralling motion 
is not possible for bubbles of radius larger than 3-6mm. ‘This agrees 
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qualitatively with the experimental observations of the gradual dis- 
appearance of the spiralling motion as the radius of the bubble increases. 
For example, Rosenberg (1950) claims that the spiralling motion ceases 
when the radius becomes equal to 2:4mm. We shall see below that there 
are indications that the figure of 3-6 mm would be reduced somewhat if 
keZy/(WQd) were not neglected. 

The ratio a/b, as given by the analysis, of the longest to the shortest axis 
of the bubble is plotted against the radius in figure +. ‘The experimental 
measurements of this ratio by Rosenberg (1950) are also shown. ‘The 
agreement is not too good, especially for the larger bubbles, but we shall 
see that the neglect of kgZy/(WQqd) leads to the value of a/b for the larger 
bubbles being somewhat too small. Another reason for the discrepancy 
may be the assumption that the bubble is exactly an oblate spheroid. 
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Figure +. [he shape of air bubbles. 


The velocity of rise as given by equations (11) and (13) is shown in 
figure 5, together with experimental results for spiralling bubbles given 
by Datta, Napier & Newitt (1950) and Haberman & Morton (1953). 
The theoretical and experimental results are in qualitative agreement but 
the theoretical values are about 20% too small. The theoretical velocity 
goes off to infinity as the radius decreases, but small bubbles do not spiral 
and the analysis would therefore not apply to them. We shall discuss below 
reasons why the calculated velocities are too small. 

We see from equation (12) that Qd/(Wy) decreases as r increases, 
being one when r is zero, about 0-75 when r = 2:3 mm, and about 0-2 when 
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r=3-6mm. Provided that the bubble radius is not too large, this result 
is consistent with the rough experimental observations that Qd/(Wy) = 1 
(see, for example, Miyagi (1925)). (It follows from (12’) below that for 
large bubbles Qd/(Wy), as predicted by (12), is too small.) 

Let us now consider the magnitude of some of the terms that were 
neglected during the course of the analysis. First, there are assumptions 
that certain geometrical quantities, such as the angle of the spiral Qd/W 
and y, are small. With typical values such as 2 = 30 radians/sec, 
d= 1-5 mm and W = 30 cm/sec, it is easily seen that all such geometrical 
quantities that have been neglected could possibly account for the dis- 
crepancies between theory and experiment, but this does not seem likely, 
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Figure 5. The velocity of rise of spiralling bubbles 


Of more importance is the neglect of kgZy/(WQd) = A, say. (The 
physical meaning of A is not altogether clear, although the requirement 
that it is small probably expresses the condition that variations of hydro- 
static pressure over the bubble surface are small compared with those of 
the hydrodynamic pressure.) When r is less than 1:5 mm, A is less than 


(-25, and its neglect is reasonable. However, if we determine y by 
equation (12), A increases rapidly with r and is of the order of unity when 
the radius is about 3-0 mm; therefore equations (12) and (13) will become 


more and more inaccurate as the bubble becomes larger. 
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If we did not neglect A, equations (12) and (13) would be 





Qd X* 
Wy 2 (12’) 
and 
 —— (1+) ' 
mR + ep ~ 2YEX*(I-Ayz—x*) | (13’) 
where 
r* Q4d2) 12 O42) -12 
Z = (ayes Sah q+ Sah. 
bs ay, ae q ty 


We cannot, of course, solve (13’) for €; but it can be shown that since 
X* > X, equation (13) gives values of ¢ which are too large, that is, the 
actual values of ¢ should be less than those shown in figure 3. It follows 
that the values of a/b should be larger than those shown in figure 4. 

We can obtain a solution of (12’) and (13’) if we suppose that A is nearly 
equal to unity. ‘This supposition is an extra relationship between the 
variables that might never be satisfied in practice, i.e. it might be inconsistent 
with equations derived from a study of other regions of the flow; conse- 
quently results obtained by its use may be inaccurate, but they should give 
an idea of the effect of not neglecting A. When A = 1, 

Q4d?/(g?y?) = O7k2Z?2/W? < |, 
and so X* = A/(1—A). From (13’), it follows that X* must be less than Z 
if r is to be real, and so X must be small when (1—A) is small. Thus, it is 
a good approximation to take ¢ as the root of X = 0, 1e. €=0-3. On 
combining equations (11), (12’) and (13’) with the condition that A is 
nearly equal to unity, we obtain X* = 2YZ(1+2Y)"7 = 0-65, r= 1-3mm 
and W = 25cm/sec. 

It can be shown further that A cannot be larger than unity, and so we 
have the result that the spiralling motion ceases when r = 1-8 mm, 
W = 25 cm/sec, a/b = 3-5 and Qd/(Wy) = 0-54. The velocity of rise is 
in closer agreement with the experimental results than the value obtained 
previously, but a/b is too large and the maximum radius for spiralling 
bubbles is too small; thus it seems likely that the supposition of A being 
nearly equal to unity is inconsistent. However, the above result indicates 
that, if A were not neglected, the effect would be to give a maximum radius 
smaller than that obtained previously, a slightly increased W and an 
increased (Q2d/(Wy) for bubbles whose radius is nearly equal to the maximum 
radius. 

Although, as we have seen, the neglect of A may account for some of 
the discrepancies between theory and experiment, it will not account for 
the discrepancies between the calculated and observed velocities of rise, 
especially when the bubbles are small. A possible explanation lies in the 
assumption that the bubble is exactly an oblate spheroid. If the bubble 
were spherical, of radius R, equation (11) would become W = 2(gR)!?; 
and when, for instance, R= 1-0mm, W would be 6:7 cm/sec. It is 
therefore clear that the calculated velocity of rise is sensitive to the assumed 
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shape of the bubble. Now the surface tension enters into the above equations 
as a product with the rate of variation along the surface of the sum of the 
principal curvatures, and if this rate were larger, then so would be the 
velocity of rise. ‘There are thus reasons for expecting that if the assumed 
shape of the bubble had greater variations of curvature than an oblate 
spheroid, as may well be the case, then better agreement with experiment 
would be obtained. But if such a shape were assumed, the mathematical 
analysis would become much more difficult, if, indeed, possible at all. 


Separation of the flow 

The analysis can be moditied in order to take some account of the 
possibility that the flow around the bubble separates. It was previously 
assumed that the flow at the front of the bubble is the same as the ideal 
How around a body of the same shape. However, there is evidence of a wake 
behind a spiralling bubble (see, for example, Miyagi 1925), and this implies 
that the flow around the bubble separates. It is still reasonable to suppose 
that the fow near the front is irrotational, but the flow will be different 
from that given by expression (3). The principal difference will probably 
be in the gradient of the velocity at the stagnation point, and it is unlikely 
that equations (+) and (5), which determine the position of the front 
stagnation point, are seriously in error. Equations (7) and (8) will still 
be correct, since they arise from the first derivatives in (6) evaluated at 
the stagnation point, and are independent of the velocity. 

However, equation (11) needs modifying since it involves ¢*V*/cs- 
(where €/cs denotes ditferentiation along the surface), and this will attect 
the coefficient of H*. It should be replaced by 

| KZ? 4¢7 > 

ae 7 Tee eee tt) 
where A is the value of 0?V?/cs* for the ideal flow divided by the value tor 
the separation flow. ‘l'ypical values of AK are probably obtained by comparing 
the ideal flow around a sphere with the ideal flow around a variety of Rankine 


(11°) 


solids, it being supposed that the separation flow around the sphere 
corresponds to the ideal flow around the solid. Values of K between one 
and two are then obtained, depending upon the shape of the solid; and it 
seems plausible to suppose that the values of A for the flow around a bubble 
will also lie in this range. 

Equation (12) will then be replaced by 


Qd ~ 
iT. = ED (12") 
Wy KZ 
and equation (13) by 
47 1+? , 
_—_— : = (13”) 


pe(+@ ~— 2YQZ—X/K) 

If A is less than 1-5, the solutions are similar to those for K = 1, which 
have already been given. The limit r = 0 still corresponds to a finite 
value of ¢, although this value is larger than 1-52 and tends to infinity as 
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A tends to 1-5. Also, IV tends to infinity as r tends to zero, but is less or 
greater than the corresponding values for K 1, according as r is less 
or greater than some value in the neighbourhood of 1-4 mm (the exact 
value of r at the ‘cross-over’ depending on A), but the values of W are 
still close to those for K E: 

However, if A is greater than 1-5, ¢ tends to infinity as r tends to zero, 
while W’ is zero when r = O and increases monotonically with 7, crossing 
the curve corresponding to A 1 in the neighbourhood of 16mm. In 
ill cases, however, there is a maximum value for the radius of a spiralling 
bubble which ts close to the maximum value that was obtained with A = 1. 
rhe solutions of equations (11’), (12”) and (13”), with a typical value 
K 1-75, are shown in figures 3 and 5. 

Of special interest is the result that the curves of W against r, when K 
is greater than 1-5, are similar in shape to the experimental curve for air 
bubbles in tap water or water containing surface-active agents obtained by 
Haberman & Morton (see figure 5), the theoretical values of Wwith A = 1-75 
being about 10°,, too small. (Closer agreement could be obtained here by 
choosing K smaller than 1-75, but this is of little significance.) “hus, the 
theory can be made consistent with the experimental observations of bubbles 
in tap water, provided that we postulate increased separation of the flow 
ind the possibility that A may be greater than 1-5. 

It has been suggested in the literature that the effect of impurities or 
surface-active substances is to make the surface of the bubble more ‘ solid’, 
thus increasing the drag and reducing the velocity of rise. ‘This suggestion 
is consistent with the above result, since if the surface of the bubble is more 
‘solid’, then, presumably, separation will occur earlier and A will be 
increased. The fact that the drag in greater also suggests that separation 
occurs earlier, since the skin friction is probably small compared with the 
form drag (at a pure air-water interface, the skin friction is zere), and an 
earlier separation of the flow usually accompanies a marked increase in 
form drag. ‘The increase in drag has been associated with the existence of 
circulation of the air inside the bubble, but it is difficult to see how 
a motion of the air can appreciably affect the motion of the water, since 
the density and viscosity of air are so small compared with those of 
vater. 

It is to be emphasized, however, that there is no direct evidence which 
relates the presence of impurities to the separation and thence to the value 
of K; and until such evidence is forthcoming, these ideas are to be treated 
with reserve. 

It is perhaps worth mentioning here that Davies & ‘Taylor (1950) 
found that the flow near the nose of a large spherical-cap bubble was close 
to the ideal flow past the sphere of which the bubble was the cap. Since 
photographs of these bubbles indicated that a bubble and its wake together 
had approximately the shape of a sphere, so that the dividing streamlines 
lay approximately on the surface of a sphere, this result is not in contradiction 
with the above concept of separation of the flow. 
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The velocities of spiralling bubbles in liquids other than water have 
been measured by Haberman & Morton (1953), and the curves of velocity 
of rise against radius are similar in shape to the experimental curve (a) shown 
in figure 5, except that the peaks are not quite so sharp. It follows from 
equations (11) and (13) that the dependence of W(T/p)-!4 on 7(T/p)-"? is 
the same for all liquids, and the data given by Haberman & Morton satisty 


this criterion with reasonable accuracy. 


4+. ‘VHE STABILITY OF A RISING BUBBLI 


The method of $3 will be applied in this section to the motion of a 
zig-zagging spheroidal bubble. (As mentioned previously, we shall obtain 
from the analysis an equation which determines the stability of the récti- 
linear motion of a spheroidal bubble.) It is necessary to take into account 
the etfects of variations with time of the shape, orientation and velocity of 
the bubble, since the motion of a zig-zagging bubble is not steady. ‘lo make 
the analysis tractable, such variations will be treated as small perturbations 
of the rectilinear motion. However, to account simultaneously for 
variations of these three quantities would still make the analysis unduly 
complicated; and, in order to find whether any simplifying assumptions 
are permissible, we shall first treat the case in which the unperturbed shape 
is taken to be a sphere instead of a spheroid. The results for a spherical 
bubble will enable us, in fact, to make plausible assumptions about the 
motion of a spheroidal bubble and to simplify the analysis considerably. 


The motion of a spherical bubble 

We calculate the velocity and pressure of the water at the front of the 
bubble using the same assumptions as in § 3, that is, supposing the flow to 
be the ideal flow around a body of the same shape as the bubble. The 
unperturbed shape is taken to be the sphere r = R, and the perturbed shape 
to be? R(1 + €)So + €,S, + €oS,), where €9Sq = €9(t), €,S, = Ox+ oy + ds, 
and €,S, = Ax?+ wy? + vs? + 2ays+28sx+2yxy, where A+ u+v=0. It is 
assumed that sufficient generality is obtained by taking the perturbed shape 
as an ellipsoid and neglecting higher harmonics. ‘The unperturbed velocity 
of rise is denoted by W = (W, 0, 0), where the x-axis is taken vertically 
upwards. Perturbations of the velocity of the centre of the bubble can 
easily be shown to be equivalent to the term ¢,S,, that is, a perturbation 
(u, v, w) of the velocity may be represented by €,S, with R*(0, 4, hb) =(u, 7, w) 

We denote the velocity potential by ©. Then ® satisfies V?® = 0, and 








--Oasr-» «x. At the surface of the bubble, c@/cn = cr/et + W.n, where 
n denotes the unit normal to the surface. It follows, after some algebra, that 
WxR® 2 RS, € 20-1 
eee ee iar a ae = 
2r’ ‘6 n+1 = 
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where the dot denotes differentiation with respect to time, and the term 
with m in the denominator is to be omitted when m = 0. The pressure is 
given by 

p (VT Dy2 1D — r , z 

~ +](V0)?+d®0 dt+gx = constant, 

e 
where d/dt = (¢/ct—We/ex). Also, the sum of the principal curvatures 
of the surface is (see Lamb 1932, § 275) 

l l l : 


— + — = > (n—1)(n+2)e,S,/R. 


a —_ = 


Hence, the value of p+7(1/R,+1/R,) at the surface of the bubble ts 
known. As in §3, we suppose that p+ 7(1 R,+1 R,), when calculated 
in this manner, is such that its first and second derivatives (differentiating 
along the surface) vanish at the point y = = = 0. 

After some algebra, we obtain the following equations: 


: 15 Wx ee 3g ail 4 
oo oe oe + Re = ve 


and the same equation with « replaced by n—v; also 
93Wy ie =| 30, 27 Wi 


20R a ton) * ae * 3 (15 








and the same equation with y and wv replaced by § and w respectivel, ; and 


9W2 g 3We, ge lun O9Wu 
$e RTDR R*IR*ER 
_ 81... ./12T 8eR\ 
+ Rey + 0 WR. 4a(— ~ 7 = 0. (16 


If the first and second derivatives of p+7(1/R,+1/R,) evaluated at 
some other point of the surface had been put equal to zero, then the same 
equations would have been obtained provided the distance between this 


point and y = s = O was of the order of «. For extra terms of this order 
would arise only if the third derivatives of p+ 7(1/R,+1/R,) evaluated at 
= = 0 were not small, whereas, since 


a l DS Wee... ; 
p+T( x * z) = <= GP pe Ls 3“) + U(e), 


the third derivatives are in fact of the order of e. 

It follows from equation (14) that the bubble shape is always stable 
to perturbations of « and ~—v, and that such perturbations are heavily 
damped. It can be shown that the deviation from axial symmetry about 
the vertical of the shape of the bubble is dependent upon these quantities. 
Equations (15) and (16) involve more than one variable, and no informtaion 
about stability can be obtained without the aid of other relations between 
the variables; how such relations are to be obtained is not clear. Of 
importance, however, is the result that y and v occur together in an equation 


S2 
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which does not involve any of the other variables. ‘The meaning of y is 
clear when the perturbation is symmetrical about the plane x, y, and eS, 
can be taken as Ax? + py? + v2*+2yxy. Then y is proportional to the angle 
between the principal axes of the distorted sphere and the coordinate axes, 
so that equation (15) involves only the sideways velocity v and the orientation 
of the bubble. 


The sig-sagging sphe roidal bubble 

As pointed out at the beginning of § 3, the assumption that a zig-zagging 
bubble is spherical is a bad one. For example, equation (16) predicts that 
W = ?(gR)'*, which gives values of W that are too small by a factor of 
three. We shall therefore assume that the bubble is an oblate spheroid. 
Further, we shall neglect changes of shape and vertical velocity, and 
consider variations of the sideways velocity and the orientation of the bubble 
only. ‘here are two reasons for this. Firstly, it seems plausible to assume, 
by analogy with the equations obtained above for a spherical bubble, that 
departures from axial symmetry are heavily damped, and that the equation 
involving the sideways velocity and orientation of the bubble is independent 
of perturbations of the vertical velocity and shape. Secondly, it is 
variations of the sideways velocity and orientation only that are directly 
observed. 

We take axes (x, y, 3) fixed in the bubble, x being the axis of symmetry, 
and x, v the vertical plane containing the zig-zag. Let y (measured in the 
positive sense from x) denote the angle between the x-axis and the vertical. 
We suppose that the bubble moves with constant vertical velocity HW’ and 
variable horizontal velocity v in the plane x, y._ It is supposed that y and 7 
are small, and their squares and products are neglected. 

Referred to the axes (x, v, s) the bubble has velocity (W, Wy +7, 0) 
and angular velocity (0, 0, — x). We use oblate spheroidal coordinates ¢, 
yo and w (see $3). Then the velocity potential ® is given by 


‘p RW1C,(1+ G)-1 — cot fy} (1 — f cot 12) + k(v + Wy) x 


i 





249 y 
a , 9\1/¢5 -9\1/9 S , 
x {—— cot!{y> (1—p?)'2(1 + C7)? « = — —- cot >cosw + 
Co( Cn + 1) ee l 
+ Ap(1 — po?) 21 + 2?) 230 cot fC -3 + (14+ 2) eos, 
where the surface of the bubble corresponds to € = , and 


A E (22+ 1)!°[30cot-Z —3 +14 2) | = BY(14+ 2722, 
Ihe pressure is given by 

pip +i(VO)?—U-V0+60/et+gH = const., (17) 
where U = (W, Wy +2, 0)+ (yy, — xx, 0) and H = xcosy+ysin x. 


To determine v, y and W, we shall use the condition that the first and 
second derivatives along the surface of p+ 7(1/R,+1/R,) vanish at the 
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point where the dividing streamline meets the surface,* this point being 
given by V® = U. 

After some algebra, it is found that the coordinates of this point are 


yA 


given by 


o= 0 sng = Wa eye 2h t Wy)- Xky), (18) 
where « = cos@, and the suffix of fj has been dropped since all expressions 
are evaluated on € = (. We put q = V®—U, and rewrite equation (17) 
as 

D ‘ ep 
: T $q° _ 4U? 1 ey, Tt gH const. 


On neglecting squares and products of y and v, we obtain 


Ww*(1— we) + ‘ga 2W(1 + ¢7)8 21 — p?)P? 


(+O (eLAyZ 
( 22 | Wy) Vp : vR[(X ayG 217) c] ), 
U? =.W?+2Why(1 + 22)"? sin 9 cosw, 


— = —k(t+Wy)LYZ(1+ 2)? coswsin 64 


+ k*y(1 + €)!*(X — 0) cos @sin 4 cos w. 


1 (1+ 2) ] l 
Rn, ky k (C24 p2)82 7 (2419+ p22 * 
We now apply the boundary condition on the pressure. ‘The first 
derivative with respect to w vanishes identically. ‘The vanishing of the 
first derivative with respect to 6 gives, after some algebra, * 


f 4T¢ 
ee 1 (2)-12 gj 
( 1)? sc) ( Co) in@ 


Si, (ng aac ee 
-Wx-(@+ Wi) irn +gytky(X-0). (19) 





‘The second derivative with respect to @ and w vanishes identically. ‘The 
second derivative with respect to @ gives 

Ww? Zz 47T¢ 

gk al 1+ 
and the second derivative with respect to w is a combination of equations (18), 
(19) and (11). 

On eliminating 6 between equations (18) and (19), we obtain 

CYZ xX 2YW? vl(YZ 2WYv 
ie” 2p * ee A ie oe 


eS: (11 
pgk*(1 + ¢?)? ) , 


k(Ly —X)+ Wd 14 ait 


(20) 

It is really only a matter of convenience which point is taken for the vanishing 

of the derivatives, since it can easily be shown that, provided the distance between 

the point chosen and the axis of symmetry is of order y, the same equations will 
result from the analysis. 
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It ¢ -- ow and k — Win such a way that k{ — R (this is the limiting process 
by which the spheroid becomes a sphere), it can be shown that equation (20) 
reduces to equation (15). 

Now equation (20) is an equation in the two variables wv and y, and it is 
not possible to proceed further without another relation between these 
variables. It is by no means clear how such a relation is to be obtained, 
although, presumably, one could be obtained from a consideration of the 
forces and couples on the bubble. However, these forces and couples 
cannot be calculated accurately without a knowledge of the motion in the 
wake. The simplest hypothesis that we can make in order to find another 
relation between zv and yx is to suppose that the actual couple on the bubble 
is proportional to the couple on a body of the same shape as the bubble 
moving with the same velocity in an ideal liquid. ‘This couple is easily 
calculated (see, for example, Lamb 1932, $ 124). Now, the total couple 
on the bubble is zero, since the mass of the air inside the bubble is negligible, 
and we hence obtain the result 

RRy + W(e+ Wx) = 9, (21) 
where R = $(2(14 ¢?)-{2Y(3 + 4¢?) — 3(1 + 2¢)}1. 

Now, the slowest natural frequency of oscillation of a stationary air bubble 
of radius R in water is (1/27)(12T/pR*)!?.. When R is about one millimetre, 
this frequency is about 150 per second; and, further, equations (14), (15) 
and (16) indicate that when air bubbles of about this size are rising freely, 
this natural oscillation is not sensibly affected and remains stable. The 
frequency of the zig-zag motion is about seven per second, which is much 
slower, and we shall confine ourselves here to perturbations which have 
a time scale comparable with the time scale + = 1/7 second, say, of the 
zig-zag motion. ‘That is, we shall neglect the natural oscillations of the 
bubble because they are too rapid and are apparently stable. ‘Then 


RR2¥/(W2y) = O(k2|(W2r2)) = O(10-3) 


(this equation represents the fact that the velocity perturbations due to the 
oscillation of the bubble are small compared with the velocity of rise), and 
(21) reduces to 

v+Wy = 0, (22) 
which is the condition that the velocity of the bubble is parallel to the axis 
of symmetry. ‘This result is similar to the assumption made in $3 that 
the axis of symmetry and the velocity of a spiralling bubble are nearly 
parallel; but, in the case of a zig-zagging bubble, there is no experimental 
evidence either to contradict or confirm (22), since the photographs of 
zig-zagging bubbles that have so far been taken are not good enough for 
this purpose. We shall see, however, that the use of (22) leads to results 


in reasonable agreement with experiment. 
On substituting (22) in equation (20) and neglecting terms which are 
of order k/(Wr) = O(1/30), we obtain 
Wy(Z—X)—gZy = 0. (23) 
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Equation (23), which is the final result of our analysis of the zig-zag motion, 
is essentially a stability equation and its implications will be discussed in 
the next sub-section. 


The stability of a rising bubble 

Equation (23) implies that x is damped exponentially when ¢ > 1-52, 
and increases exponentially when ¢ < 1-52, since Z—X is a monotonically 
decreasing function of ¢ which is zero when ¢ 1:52. For ¢ nearly equal 
to 1-52, equation (23) gives a large value to y/y and is inconsistent with the 
approximations by which it is derived. Hence, tor ¢ nearly equal to 1-52, 
equation (22) does not hold, or it is not valid to neglect changes of shape and 
vertical velocity, or, as seems most likely, a linearized theory is not valid 
and second-order terms should be taken into account. ‘The exponential 
increase of y when ¢ < 1-52 also implies that a linearized theory is not valid 
for long in this case. 

Thus, the analysis of the motion near the front of the bubble does not 
give us any information about the zig-zag motion, but gives the result that 
if the bubble is more oblate than the spheroid corresponding to € = 1-52, 
then it is unstable to small disturbances and will move in an irregular manner, 
although the above equations do not indicate in what way it will move. 
(It will be remembered that in §3 we found ¢ < 1-52 to be a necessary 
condition for the spiralling motion to occur.) In other words, we have found 
the value of ¢ for which the rectilinear motion of a spheroidal bubble becomes 
unstable. It is true that we have considered, in effect, only disturbances of 
sideways velocity and orientation which are symmetrical about a plane, but 
this is probably sufficient, both for the reasons given before in the sub-section 
on the zig-zagging spheroidal bubble and because the experimental observa- 
tions of § 2 show that the effect of the first instability is to produce a zig-zag 
motion ina plane. ‘lhe possible influence of the motion in the wake on the 
occurrence of the instability and the mechanism of the zig-zag motion will 
be discussed later. 

Let us now compare the above prediction with the experimental data 
for the critical bubble, i.e. the smallest bubble that is unstable. From 


equation (11) with ¢ = 1:52, T = 74 dynes/cm and g = 981 cm/sec”, we 
obtain, as the relation between the velocity and radius of the critical bubble, 

W?2 = 223k! + 81k, (24) 
where r = 1:7k (using r = k{((1 + @*)}!3), k is measured in millimetres and 


W incm/sec. ‘The observed values are r = 0-7 mm and W = 30 cm/sec 
approximately. On putting r = 0-7 mm in (24), we obtain W = 24 cm/sec, 
(W = 30 cm/sec corresponds to r = 0-4mm). This velocity is somewhat 
lower than the experimental value, as was also the case in $3 for the 
calculated velocity of a spiralling bubble, and the reasons that were put 
forward there to account for the disagreement will apply here with equal 
force. It is also worth mentioning that when ¢ = 1-52, the ratio of the 
major to the minor axis of the bubble is 1-2, which agrees well with the 
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experimental value for a bubble of radius 0-7 mm given by Rosenberg 
(1950) and shown in figure 4. 

Further, on neglecting pgk?(1 + (*)? (47), which is usually O(10 ') for 
the critical bubble (this is equivalent to neglecting the second term on the 
right-hand side of (24)), equation (11) gives the value of the critical Weber 
number 2Hrp 7 at which the instability appears. We have 

2Wrp 872243 ? 
[= e) = qe 1-03, 
/ cr 52 
whereas the experimental value for air bubbles in water ts about 1-7. 

For the sake of completeness, we shall also investigate briefly the ettect 
of separation of the flow on the relation between the velocity and radius 
of the critical bubble. As before, only equation (11) requires modification, 


and on replacing it by (11’) and repeating the analysis, we obtain 
Wy(Z—-NX K)—gZy = 0. (23’) 


The motion becomes unstable when € is equal to the root of Z— XK = 0, 
and trom (11'), we obtain the required relation between the velocity and 
radius of the critical bubble. It is found that for a given radius, the velocity 
decreases as A increases. For example, W 14cm sec corresponds to 
y= 0-7 mm when A 43. ‘These figures are consistent with the 
observations of $2 that bubbles in tap water start to zig-zag when their 
radii are about 0-S mm, even though their velocities of rise are much 
smaller than bubbles of the same size in filtered or distilled water. However, 
if KA > 3.2, Z—NX K is always positive, and this implies that if the motion 
is such that A 3,2, then it is unstable for all ¢. 


The mechanism of the zig-zag motion 

We have not so far considered the motion in the wake of the bubble, 
since no way could be found of overcoming the mathematical difficulties. 
However, there is some experimental evidence on the wake, which we 
shall now describe, and it will be seen that it suggests that the zig-zag motion 
and the structure of the wake are in some way related. 

Figure 6 shows the dependence of the drag coefficient Cp = (8/3)(gr/ W*) 
upon the Reynolds number 217». ‘This curve was constructed from data 
on the velocity of rise of zig-zagging bubbles given by Datta, Napier & 
Newitt (1950). (Haberman & Morton (1953) give a curve for spiralling 
bubbles which is very similar.) It will be noticed that the transition to 
zig-zag motion is accompanied by a rapid increase in the drag coefficient. 
All the drag is form drag, there being no skin friction since the boundary 
condition at an air water interface is zero tangential stress instead of zero 
velocity of slip, and an alteration in the behaviour of the drag coefficient 
must be associated with a marked change in the structure of the wake. 
‘There is no direct evidence on the wake behind bubbles, but there is some 
on the wake behind solid three-dimensional bodies, and this we shall now 


describe briefly. 
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Marshall & Stanton (1930) observed the flow past a circular disc which 
was normal to the direction of flow, and found that the wake was steady if 
the Reynolds number was less than 200. When the Reynolds number 
became equal to 200, a periodic discharge of vorticity, symmetrical about 
a plane perpendicular to the disc, was observed. Moeller (1938) observed 
the wake behind a solid sphere over a large range of Reynolds numbers. 
He found that for Reynolds numbers less than 170, the wake was steady 
and laminar, but for Reynolds numbers greater than 200, the wake was 
no longer steady but oscillated periodically. lhe amplitude of the 
oscillation increased with the Reynolds number, being initially quite small, 
until the wake became turbulent at a Reynolds number of about 1000, 
the turbulent wake still showing signs of a periodic structure. ‘The 
oscillation of the laminar wake was symmetrical about a plane and there was 
no evidence at all of a spiral wake. 


RECTILINEAR MOTION ZIG-ZAG MOTION 
a a aa 


DRAG COEFFICIENT 
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Figure 6. Dependence of drag coefficient on Reynolds number for zig-zag motion 


Moeller gives the frequency of oscillation of the wake N for Reynolds 
numbers larger than 800; for instance, when 2Wr/v = 800, Nr/W = 0-12. 
Marshall & Stanton give Nr/W = 0-06 for a Reynolds number of 
200 only. 

The value of Nr/W for a zig-zagging bubble is about 0-02 at Reynolds 
numbers which are greater than 400, this being the Reynolds number of 
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the critical bubble. However, we should expect the wake behind a bubble 
to be somewhat ‘ weaker’ than the wake behind a solid body, owing to the 
different boundary conditions and consequent reduction in the rate of 
vorticity production; and it seems likely that an oscillation, or instability, 
of the wake of a bubble would first appear at a higher Reynolds number 
and have a smaller frequency than the oscillations of the wake of a solid 
body. It is therefore possible that the zig-zag motion is related to an 
oscillation of the wake or a periodic discharge of vorticity trom behind the 
bubble, in so far as these may be distinguished in the light of present 
knowledge. 

We therefore suggest that the observed zig-zag motion is caused by the 
interaction of an oscillating wake and the instability of the motion near the 
front of the bubble that occurs when the bubble is sufficiently oblate 
(i.e. when ¢ is less than 1-52). It is unlikely that the wake oscillation by 
itself could produce the zig-zag, since similar motion has not apparently 
been observed with solid bodies. A familiar phenomenon that might seem 
to be relevant is the motion of a falling leaf; however, this has a different 
character, since the leaf glides through the air with its velocity in the plane 
of the leaf, whereas a bubble moves with its velocity perpendicular to its 
equatorial plane. Also, the argument that the forces arising from the 
wake oscillation are not sufficient to move a solid body sideways but would 
be sufficient to move a bubble is not valid, since a bubble has, in fact, a 
virtual mass, this being the mass of the water that has to be accelerated 
when the bubble accelerates. It may be that the instability at the front 
of the bubble occurs first and ‘triggers’ off an oscillation of the wake. 
The marked change in the behaviour of the drag coefficient makes this 
seem likely, but the actual details of how the zig-zag motion is produced 
are not clear. 

This suggested mechanism also fits the data and analysis for bubbles 
in tap water, unless A is greater than 3 and the front of the bubble is unstable 
for all ¢. In this case, we should expect the motion to be irregular until 
the Reynolds number is high enough for the oscillation of the wake to be 
possible. (It will be remembered that the analysis is only valid when the 
Reynolds number is sufficiently large for the assumption of inviscid flow 
near the front of the bubble to hold.) There is no evidence of this 
happening, and A is presumably less than 3. This does not necessarily 
contradict the supposition in §3 that K for spiralling bubbles may be 
greater than 3, since A will depend on the Reynolds number and the nature 
of the motion. 

Finally, there is no evidence concerning the causes of the spiralling 
motion, except that the experimental observations make one wonder 
whether they have anything to do with the onset of turbulence in the wake. 


The author wishes to thank Dr G. K. Batchelor for suggesting this 
problem to him and for much helpful advice and encouragement, and the 
Department of Scientific and Industrial Research for a maintenance award. 
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SUMMARY 

Results of an experimental investigation of the heat loss of tine 
heated wires immersed in a supersonic stream at right angles to the 
How direction are presented. ‘The measurements show that the 
heat loss of the wire is independent of the free-stream Mach number 
for values of the latter between 1-3 and4-5. Since the wire is always 
in the wake of a detached shock wave, the streamlines in the 
neighbourhood of the wire pass through a normal shock wave. 
‘he Reynolds number Re,, based on conditions behind a normal 
shock, thus becomes the characteristic parameter for the heat 
transfer. ‘lhe measurements covering a Reynolds number range 
of 3 to 220 show the existence of two flow regimes. For Re, > 20 
the Nusselt number is a linear function of the square root of the 
Reynolds number, and the equilibrium temperature is nearly 
independent of Re,. For Re, — 20 the Nusselt number decreases 
more slowly, and the equilibrium temperature rises sharply with 
decreasing Reynolds number. 


INTRODUCTION 


‘The determination of the drag and heat transfer of a cylinder in a viscous 
compressible flow is a classical problem in fluid mechanics. At present, 
no theoretical solution valid for a wide range of Mach and Reynolds numbers 
exists. Although several accounts of experimental work on the heat 
transfer problem have been published (Kovasznay 1950; Lowell 1950; 
Spangenberg 1955; Stalder 1952; Stine 1954; Weske 1943), they disagree 
without exception in the manner of presenting the results; and, in cases 
where sufficient information is avatlable for direct comparison, the results 
also disagree in absolute values. It was felt advisable, therefore, to carry 
out a series of experiments that might not only resolve the discrepancies 
found in existing measurements, but, more importantly, might throw some 
light on basic questions concerning (a) the principal parameters involved, 
(6) whether or not heat transfer from blunt bodies in compressible flows is 
basically different from that in low-speed flows. 


With regard to (a), the point of view of the investigation was influenced 
by some recent experiments on the flow field around blunt bodies at super- 
sonic speeds (Stine 1954; Walter 1953). These experiments have indicated 
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that there is a tendency for the How field behind the detached shock wave 
in the vicinity of the blunt body to approach a fixed pattern as the free-stream 
Mach number is increased. If the concept of the frozen flow field in the 
neighbourhood of the body is indeed correct, it would greatly simplify 
the heat transfer problem. It would mean that the local quantities 
characteristic of the frozen pattern, rather than the free-stream parameters 
(M,, Re..), the 
the local parameters should be based on conditions behind the detached 
shock wave (Re, .W/,). 
numbers, .W, can be considered constant, the Reynolds number Re, can 
[t is shown 


control heat transfer. It is therefore suggested that 


Since, for sufficiently high free-stream Mach 


thus be expected to be the principal parameter of the problem. 
in the paper that this is in fact the case. 

With regard to (4), at present little is known about the flow field between 
a blunt body and the shock wave in front of it in flows at high Mach number 
and low Reynolds number. Although the present experiments were not 
designed to investigate this problem, it was possible to obtain some informa- 
tion relevant to it. ‘Ihe measurements provide a lower Reynolds number 
limit (Re, ~ 20) above which the usual boundary layer approximations 
are expected to be valid. At lower Reynolds numbers the measurements 
indicate that the flow field changes, but no definite statement can be made 
as to the nature of this change, except that here too the parameter Re, 
correlates the results satisfactorily. 

The experiments were carried out under joint sponsorship of the 
Department of the Army, Ordnance Corps(under Contract No. DA-04-495- 
Ord18), and the Department of the Air Force. ‘The authors wish to 
acknowledge valuable discussions with Drs Peter P. Wegener, Leslie Mack 


and H. W. Liepmann. 


SYMBOLS 


a speed of sound, cm/sec. 


Re pud , Reynolds number. 
¢, specific heat at constant T 


temperature, K. 





pressure, calgram™!?° C?. T; local stagnation temperature in 
d wire diameter, cm. boundary layer. 
i wire current, amps. t temperature, © C. 
h heat-transfer coefficient, u_ velocity, cm/sec. 

calenr* sec ** C™, v distance from flat plate leading 
k thermal conductivity of air, edge, cm. 

calcar*sec*” C—*. y distance from flat plate surface, 
1 wire length, cm. cm. 
M_ u/a, Mach number. x temperature coefficient of 
Nu hd/k, Nusselt number. resistance, © C~!. 
p pressure, cm Hg. y _ ratio of specific heats. 
p, Pitot pressure, cm Hg. n (y/x)\/(Re,) boundary layer 
Pr c,p/k, Prandtl number. parameter. 
yg heat loss from wire, cal/sec. pe absolute viscosity, poise. 


R_ wire resistance, ohms. p 


density, g/cm?. 
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7 temperature loading, equilibrium or unheated wire 
(T,,—T,)/T.. condition. 
stagnation condition. 
Subscripts heated wire condition. 
GO aoc. undisturbed or  free-stream 
2 conditions behind detached condition. 
normal shock wave. 


EXPERIMENTAL EQUIPMENT 

Wind tunnel 

All measurements were performed in the 20in. supersonic wind tunnel 
of the Jet Propulsion Laboratory, California Institute of Technology. 
This tunnel is continuously operated with dry air by suitably staged 
centrifugal compressors, the supply pressure being variable over a range 
():2 to 4:3 atmospheres and increasing with Mach number. ‘The supply 
temperature for most experiments was between 30 and 40° C. ‘The tunnel 
supply temperature is known to be uniform across the wind tunnel settling 
chamber within 1° C. Mach numbers in small increments between 
1-3 and 5-0 can be obtained in the test section with a servo-driven stainless 


steel flexible-plate nozzle. 


Hot wires 

The wire material used was an alloy of 90°, platinum and 10°, rhodium. 
The nominal wire diameters were 0-00127 cm (0-00051n.) and 0-00038 cm 
(0-00015 in.) as specified by the manufacturer. No attempt was made to 
determine wire diameters independently since use of the nominal values 
yielded consistent data in the final results. ‘The length to diameter ratio, 


or aspect ratio, of the wires was approximately 400 for the larger diameter 


wires and 550 for the smaller ones. ‘The wire holder shown in figure 1 
(plate 1) was made of two stainless steel wedge-shaped prongs insulated 
electrically from one another by a layer of glass cloth bonded to the metal 
with a Teflon bonding agent. ‘The wires were attached to the prongs of 


the holder with soft solder. 


Heating circuit 

The direct current for heating the wires was supplied by a 24-volt 
storage battery through a circuit including a Wheatstone bridge for the 
measurement of the wire resistance. A Brown null indicator showed 
the bridge balance. A Rubicon portable precision potentiometer was used 
to measure the potential drop across a precision resistor of 1 ohm in series 
with the hot wire. 


Pitot probe 

In order to determine the local Mach number accurately, a Pitot probe 
was placed close to the hot wire. The Pitot pressure was measured on 
a mercury micromanometer to an accuracy of +0-1mm mercury. 
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Figure 2. Traversing mechanism. 
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Traversing mechanism 

The boundary layer survey presented was made on a smooth flat plate 
which spanned the working section of the tunnel. A traversing mechanism 
(figure 2, plate 1) fitted in the plate could be positioned to an accuracy of 
+ 0-01 mm. 

[}XPERIMENTAL PROCEDURE 

For the determination of the local flow conditions, the standard methods 
of recording the tunnel stagnation pressure, stagnation temperature, 
and the local Pitot pressure were used. ‘The accuracy of the measured 
Mach numbers is believed to be within }°,, and the values of the Reynolds 
numbers are accurate within 2°,,. 

In order to obtain the heat loss coefficient, or Nusselt number, a 
measurement of wire resistance and current is required, together with a 
knowledge of the heated and unheated wire temperatures. Since measure- 
ments of resistance and current can be carried out with great precision 
by well-known methods, only the accurate establishment of the wire 
temperatures necessitated special effort. This involved the measurement 
of the thermal coefficient of resistance «, carried out in the following manner. 
The resistance of several samples of the 0-00127cm and 0-00038 cm 
diameter wires was measured at 0° C, at room temperature, and near the 
temperature of boiling water. The coefficient «, based on 0° C, was 
computed from R = R,(1+t) with R being the resistance of the wire at 
some temperature ¢ C. It was considered that within the temperature 
range 0° C to 300° C covered in the present experiment, the use of a linear 
resistance-temperature relationship was acceptable, A typical calibration 
curve is shown in figure 3. It was found that for the large diameter wire 
a = 0-00175/° C+ 3°,, and for the smaller wires « = 0-00166/° C+ 13%. 
The percentages refer to the maximum variations of the thermal coefficient 
of resistance for various samples taken from the same spool of wire. ‘The 
above values of « were used in reducing all the heat loss data except at Mach 
numbers 1-33 and 4:54. In these cases a wire from a different spool with 
a = ()-00169/° C was used. 

In principle, having the information on « and R, and a resistance measure- 
ment for a given flow condition, the wire temperature could be measured. 
Unfortunately, it was found that during several measurements the wire 
stretched permanently, presumably due to the wind tunnel starting loads, 
or even due to continuous high air loads. The wire stretching increased 
the resistance by several percent. ‘This fact was established by a resistance— 
temperature recalibration of a few wires after several hours of exposure to 
the air stream. This occasional stretching made it impossible to accurately 
infer the wire temperature from its measured resistance. It was therefore 
decided to make a number of careful separate measurements of the 
equilibrium (unheated) resistance of wires exposed to the free stream in 
the entire Mach number and Reynolds number range, and for this purpose 
some specially constructed wires were used. The length of these wires 
was somewhat smaller than that used for the later heat loss measurements 
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in order to increase the probability of test survival. (Calculations indicated 
that in the case of unheated wires, the correction due to heat loss through 
the prongs of the holder was within the experimental scatter.) Furthermore, 
the resistance of all the wires used in this set of experiments was calibrated 
against temperature before and after the measurements. If the resistance 
had changed more than a few tenths of one percent due to exposure to the 
air stream, the results were discarded. With this technique it was possible 
to obtain wire equilibrium measurements repeatable within + 1". 
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Figure 3. Temperature-resistance calibration. 
Ihe procedure for the heat loss measurements may be described as 


follows. Withthe wire inthe free stream, its ‘ cold’ resistance was measured, 
that is, the resistance of the wire when there is no heating current in it, 


as in measurements of the wire recovery temperature. Next the wire was 
heated by passing a current through it until its resistance reached some 
predetermined value, and this current was measured. ‘The sequence was 
then repeated for a different predetermined value of resistance. 
Simultaneously with these measurements, the tunnel supply pressure and 
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temperature were recorded, and the Pitot pressure in the vicinity of the 
wire was measured. The entire procedure was repeated at several tunnel 
pressure levels at a given Mach number, so as to give as wide a variation of 
Reynolds number as possible. In turn the Mach number was also varied 
ver the entire range of the tunnel. Figure 4 shows a typical set of these 
measurements at M,, = 3-05. 
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DATA REDUCTION 
The measurement of the wire heat loss involved the recording of the 
following quantities. 
(a) The flow parameters: stagnation pressure p,, Pitot pressure P,, 
and stagnation temperature 7,. From these the local Mach and Reynolds 
numbers were obtained. The flow conditions on which the Reynolds 
number was based will be discussed later. 
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(6) The wire parameters: current 7, unheated or equilibrium resistance 
R,, heated wire resistance R,. From the values of R, measured at 
various Mach and Reynolds numbers, the equilibrium wire temperatures 
could be computed using the relation 
R,—R, 
Ry 

where «and R, were obtained from the resistance-temperature calibrations 
Thus the functional relationship 

a, 

7 = f(M.,, Re.) (2) 


tT - 279°C, (1) 


could be established. 

Only wires that survived a complete run and showed no change in R, 
on recalibration were used in the determination of this function. Once 
established, this relationship permitted the computation of 7’, from the 
known values of 7,, 17, and Re, in subsequent measurements. 

The heat loss from the wire was expressed in the torm of a Nusselt 


number 


Vu 
where 
Y 
~ wld(T,—T,) 
is the heat transfer coethcient foracylinder, In terms of electrical quantities, 
y 0-23977R,.R 
“i Sy 7) * see -R.)[1+2(T,, — 273)]’ 


where ():239 is a conversion factor from cal sec to watts. ‘The temperatures 
on which the values of the heat conductivity of air were based will be 
described later. 

Since R, was directly measured and 7, was obtained from the previously 
established experimental relation (2), any change in R, due to stretching 
of the wire during the supersonic measurement could be immediately 
detected, and the correct value of R, determined. ‘This change was always 
less than 3°, of the initial value of Ry. In this simple correction it was 
assumed that changes in resistance resulted from uniform changes in the 
length and diameter of the wire. 

Due to the tact that the hot wire is of finite length and the holder is at 
a lower temperature ‘than the wire itself, there is heat conduction to the 
holder as well as heat loss to the air stream by forced convection. ‘This 
means that the measured heat loss per unit length for a wire of finite length 
is different than it would be for a wire of infinite length in the same air 
stream with the same heating current. ‘The correction for this end loss 
effect was given first by King (1914). The same method was used here 


adopting the technique of computation of Kovasznay (1950). The 
corrections were of the order of 5°,, of the measured Nusseit number. 
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RESULTS AND DISCUSSION 
The principal aim of the present investigation was, first, to find the 
important flow parameters that govern the heat loss from a heated wire 
in supersonic flow, and, second, to obtain an empirical relationship between 
these flow parameters and the measured heat loss. 
In general it can be shown from dimensional arguments that the heat 
loss, usually expressed in terms of a nondimensicnal Nusselt number, 
depends on the following quantities (see Kovasznay 1950) : 


Nu,, = Nu(Re,,, M,,, Pr... y, 7, l/d). (3) 
\ll the parameters are referred to free-stream conditions, which are 
usually given in the problem. ‘lhe Prandtl number and the ratio of the 


specific heats of the gas were very nearly constant in the present experiments. 
The aspect ratio //d enters into the problem mainly because of conduction 
losses through the wire supports. .A suitable small correction described 
earlier was made to eliminate this effect. ‘The above functional relationship 


simplifies, therefore, to 


Nu, = Nu(Re,,, M.,, 7). (4) 
‘he heat transfer measurements presented in this form exhibited a 
systematic variation with free-stream Mach number. ‘The question 


arose whether, with a more judicious choice of parameters, further simpli- 
fication could be attained. 

In this connection the experiments of Walter & Lange (1953) are of 
interest. ‘Chey measured the local recovery temperature and the pressure 
distribution around an insulated cylinder in supersonic flow. ‘Thei 
measurements indicated that the flow conditions behind the detached shock 
wave and in the vicinity of the cylinder were very nearly independent of the 
(2-5 WV =< 5-0). 


free-stream Mach number in the range investigated (2 
Although the Reynolds numbers of these experiments were much higher thati 
those met here, the cylinder boundary layer was laminar, and their conclusions 
can be expected to hold also in the lower Reynolds number range. ‘This 
immediately suggests that the parameters in (+4) should be expressed in 
terms of the conditions beyond the detached normal shock wave, rather 
than the undisturbed flow conditions. Accordingly, (4) can be rewritten 
} 
, Rg ; [Ly - 
Nu, — = Nul Re,.—, M,,7}. (5) 
ihe ~ ph 7) 
The ratios of heat conductivity and viscosity depend very nearly on the 
temperature ratio 7/7, which in turn is a function of M, only. ‘Therefore 
(5) becomes Nu, = Nu(Res, My, 7). 
As the supersonic free-stream Mach number increases, M, approaches a 
constant. ‘Thus, 
Nuy = Nu( Res, 7) for M,, > 1. (6) 
All the measurements were expressed in terms of these parameters and are 


discussed below. 
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Figure 5 shows the results of measurements of the equilibrium 
temperature attained by the unheated Wire at various Mach and Reynolds 
numbers. It is seen that for a Reynolds number Re,, larger than approxi- 
mately 20, the equilibrium temperature variation with Re, and M is 
negligible, being within the limits of the experimental accuracy (+ 1 
Unfortunately, the Reynolds number range is not very wide at the higher 
Mach numbers; but the results of Walter & Lange (1953) strongly suggest 
that the equilibrium temperature has a constant value. Below a Reynolds 
number of 20, on the other hand, a sharp rise in the equilibrium temperature 
was observed, and values above the stagnation temperature were measured. 
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Figure 5. Variation of equilibrium temperature. 


The results of the heat loss measurements are shown in figure 6. The 
Nusselt number is plotted as a function of the Reynolds number with 
the temperature loading as a parameter. Both Nusselt and Reynolds 
numbers are expressed in terms of conditions behind the normal detached 
shock wave. ‘The points on the figure were directly read from plots of the 
type shown in figure +. It is seen that the measurements presented in this 
manner do not exhibit a dependence on the free-stream Mach number in 
accordance with equation (6), even for a Mach number as low as 1:3. 
Further, provided Re, > 20, the Nusselt number varies with the square 


root of the Reynolds number—a relationship well-known for low speed 
Hows. Again, below a Reynolds number of about 20 (4/(Re.) ~ 4:5) 
the square root relationship breaks down, and the heat loss decreases at 
a slower rate with decreasing Reynolds number. 
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From these results it is apparent that the predominant parameter of the 
problem is the Reynolds number. ‘They also substantiate the conjecture 
made earlier that the Reynolds number should not be based on the 
undisturbed free-stream conditions, and that the ‘apparent free stream’ 
of the cylinder was the flow behind the detached shock wave. For con- 
venience, and because in front of the body the shock wave was in fact very 
nearly normal, the Reynolds number was based on conditions behind 
a normal shock wave. 
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[he measurements can evidently be discussed in terms of two Reynolds 
number ranges, a Reynolds number of approximately 20 being the dividing 
line. 


High Reynolds number flows: Re, > 20 
\s pointed out earlier, in this range the heat loss coefficient varied with 
the square root of the Reynolds number for a given temperature loading. 
Such a relationship was found to exist in low-speed subsonic and transonic 
streams (King 1914; Lowell 1950; Stalder 1952), and is characteristic 
of flows where boundary layer approximations are applicable. ‘The results 
iggest, therefore, that above Re ~ 20 the usual boundary layer approxi- 
ition could be used in the theoretical treatment of this problem. 
Spangenberg’s data (1955) for MW, = 1-25 to M,, = 1-9 are found to 


‘ree with the present results within the scatter of his experiments. 
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Kovasznay and Térmarck (1950) were the first to point out the independence 
of the wire heat loss and the free-stream Mach number, but their values are 
too high by about 20 to 25”,,. 


Low Reynolds number flows: Re, < 20 

In this range the measurements indicate that (a) the equilibrium 
temperature of the wire increases rapidly with decreasing Reynolds number 
(figure 5), and exceeds the stagnation temperature, and (6) the heat loss 
decreases at a slower rate with decreasing Reynolds number (figure 6). 
The former effect was first noticed by Stalder, Goodwin & Creager (1952). 
‘They explained it from the point of view of kinetic theory, and related the 
equilibrium temperature to the ratio of the mean free path A and cylinder 
diameter. It is interesting to note that the Reynolds number Re, is 
equivalent to d/A, the reciprocal of the parameter used by Stalder, Goodwin & 
Creager, provided the viscosity relation of kinetic theory is used and the 
assumption M 1 is made. 

At this point one may inquire if, in the region where the mean free path 
and cylinder diameter are of the same order of magnitude, conditions can 
still be adequately described from the continuum mechanics point of view, 
and if it is in fact necessary to formulate the problem in terms of the kinetic 
theory. At present there is not enough evidence for a definite answer. 
However, from the fact that in the present experiments one Reynolds number 
alone—without additional parameters—correlates the measurements in 
this range as well as in the high Reynolds number range, one may expect 
continuum flow theories to be valid for the flow regime in question. 


CONCLUSION 
Measurements of the heat loss of fine cylinders in supersonic flows 

within a Mach number range 1-3 to 4:5 indicate the following results : 

(1) The principal independent parameter of the problem is the Reynolds 

number based on conditions behind the detached normal shock wave. 

(2) For Re, > 20, the cylinder equilibrium temperature is independent 

of both Reynolds number and free-stream Mach number and has a value 

T,/T, = 0-95+0-01; the Nusselt number Nu, varies with the square 

root of the Reynolds number for a given temperature loading and is 

independent of the free-stream Mach number. 

(3) For Re, < 20, the equilibrium temperature rises sharply with 

decreasing Reynolds number; in fact it exceeds the stagnation tempera- 

ture at about Re, <5. The square root relationship between Nusselt 

and Reynolds numbers does not hold; in this region also, the free-stream 

Mach number is not a parameter of the problem. 


APPENDIX 
MEAN FLOW MEASUREMENTS IN A SUPERSONIC LAMINAR BOUNDARY LAYER 
BY MEANS OF A HOT WIRE 


The results described in this paper suggest the use of a hot wire for 
measurements of some mean flow quantities in an unknown supersonic 
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flow field. As an example, hot-wire measurements were made in a laminar 
boundary layer on a flat plate 2+cm from the leading edge at a stagnation 
pressure of 24-0 cm Hg and a stagnation temperature of 285° K. ‘The free- 
stream Mach number was 3-05. ‘The measurements consisted of recording 
the unheated wire resistance and wire heat loss for a constant temperature 
loading at different distances from the flat surface. It is shown that by 
these two types of measurements and by the use of the results of this paper 
in a modified form, the mass flow pu and local stagnation temperature 7,’ can 
be obtained. Furthermore, the assumption of constant static pressure across 
the boundary layer allows the calculation of all the other flow parameters. 

The results were expressed in a more convenient form the following way. 
It was shown that, for a given temperature loading, 


Nu, = Nu(Re,). 


Che tluid properties, heat conductivity and absolute viscosity are based 
yn conditions behind the normal detached shock wave which are not known 
a priori. The above relationship can be written 


: k k . L r 
Nu,— - = f (Re, a 


Since both the heat conductivity and viscosity are functions of the 
temperature only, the ratios k,,/k, and yw, 4, depend very closely on 7,,/T,, 
which is a constant for a given temperature loading. Moreover, it is found 
that the ratio 7, 7, is very nearly unity within the Mach and Reynolds 
number range of the experiments. ‘The law of heat loss from the hot wire 
can therefore be written in the form 


Nu,, = Nu(Re,,, 7). 


Figure 7 shows the measurements evaluated in this form. ‘lhe points 
on this figure were obtained by cross-plotting from curves of Nu,, versus 
+ for constant values of Re,,, similar to plots shown in figure 4. It might 
be mentioned that the heat loss law in a form 


Nu, = Nu(Re,, 7) 


could also have been used. 
Furthermore, the variation of equilibrium temperature with Reynolds 

number was put inthe form 7 
Fr = f(Re,). 

The method of data reduction may be described as follows. 
(1) From the measurement of Re the equilibrium temperature is calculated 
by means of the relation 
R,-R, 


* xR, 


+ 271°C, 
(2) From the heat loss measurements Vu,, is computed; and, using 
figure 7, the mass flow pu can be obtained, since wire diameter and the 
viscosity based on wire temperature are known. 
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(3) From the relation T ; 
— = /f(Re,) 
T. 
the local total temperature is found. 
(+) Finally, by assuming constant static pressure across the boundary 
laver, the usual determination of Mach number, mean velocity and 
density can be made. 

Figure 8 shows the final results. The solid lines are theoretical values 
obtained by using the method of Klunker & McLean (1953).* From 
the Mach number distribution it is seen that the measured points below 
a Mach number 1-3 deviate significantly from the theoretical values. Since 
the heat loss variation shown in figure 7 is not expected to hold in this range, 
the discrepancy is not surprising. ‘Therefore, results obtained near the 
wall where ./ < 1-3 are not plotted in the other distributions shown on 
this figure. The figure indicates a very good agreement between theory 
and experiment in the region considered (.W > 1-3). It is believed that 
the useful range of the hot-wire method employed here could be extended 
to lower Mach numbers using the results of Spangenberg (1955) and by 


devising a more elaborate technique. 
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On the propagation of weak shock waves 
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SUMMARY 


A method is presented for treating problems of the propagation 
and ultimate decay of the shocks produced by explosions and by 
bodies in supersonic flight. ‘The theory is restricted to weak 
shocks, but is of quite general application within that limitation. 
In the author’s earlier work on this subject (Whitham 1952), only 
problems having directional symmetry were considered; thus, 
steady supersonic flow past an axisymmetrical body was a typical 
example. ‘The present paper extends the method to problems 
lacking such symmetry. ‘lhe main step required in the extension 
is described in the introduction and the general theory is completed 
in $2; the remainder of the paper is devoted to applications of the 
theory in specific cases. 

First, in $3, the problem of the outward propagation of 
spherical shocks is reconsidered since it provides the simplest 
illustration of the ideas developed in $2. ‘Then, in $4, the theory 
is applied to a model of an unsymmetrical explosion. In $5, 
a brief outline is given of the theory developed by Rao (1956) for 
the application to a supersonic projectile moving with varying 
speed and direction. Examples of steady supersonic flow past 
unsymmetrical bodies are discussed in $$6 and 7. The first is 
the flow past a flat plate delta wing at small incidence to the stream, 
with leading edges swept inside the Mach cone ; the results 
agree with those previously found by Lighthill (1949) in his 
work on shocks in cone field problems, and this provides a valuable 
check on the theory. ‘The second application in steady supersonic 
flow is to the problem of a thin wing having a finite curved leading 
edge. It is found that in any given direction the shock from the 
leading edge ultimately decays exactly as for the bow shock on 
a body of revolution; the equivalent body of revolution for any 
direction is determined in terms of the thickness distribution of 
the wing and varies with the direction chosen. Finally in §8, 
the wave drag on the wing is calculated from the rate of dissipation 
of energy by the shocks. ‘The drag is found to be the mean of the 
drags on the equivalent bodies of revolution for the different 
directions. 


* Now at 
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1. INTRODUCTION 

In a previous paper (Whitham 1952), a method was developed for 
determining the symmetrical shocks which are produced, for example, 
by an axisymmetric body in steady supersonic flow or by a spherically 
symmetric explosion. Due to the directional symmetry in these problems, 
the flow quantities depend only on two independent variables. In this 
paper, the method is extended to deal with problems which involve more 
than two independent variables. For shocks produced by bodies in super- 
sonic flight, such problems arise when the body is unsymmetrical or when 
the velocity is not uniform; they also arise in explosion problems when 
the initial shape of charge is not spherical. 

The theory is restricted to weak shocks and for that reason the applications 
to explosion problems will be of practical value only at distances from the 
explosion which are sufficiently large for the shocks to be weak. However, 
in formulating the basic ideas of the theory, it is convenient to start with 
the problem of a weak explosion and, in fact, to consider the following 
simplified model. A region of arbitrary shape, bounded by a surface S, 
contains gas at a pressure higher than that of its surroundings; initially 
the gas is at rest with uniform pressure but at time ft = 0 it is suddenly 
released. According to the theory of sound, the wavefront carrying the 
first disturbance outwards from the explosion moves along the normals 
to the surface S with velocity a), where a, is the constant sound speed in the 
undisturbed gas surrounding the explosion. ‘These normals are the 
orthogonal trajectories of the successive positions of the wavefront and are 
known as ‘rays’; in a sense, the rays are the carriers of the disturbance. 
Moreover, the appropriate solution to this problem in the theory of sound 
predicts the magnitude of the disturbance, and in particular the variation 
in the magnitude of the pressure jump at the wavefront as it moves out 
along a ray. However, near the head of the wave, the law of variation of 
the amplitude of the disturbance takes a simple form which can be deduced 
quite generally from the approximation of ‘ geometrical acoustics’, without 
appeal to the detailed solution in the full theory. For, in certain circum- 
stances, the energy propagated down a narrow ray tube formed by a bundle 
of neighbouring rays is conserved; that is, reflection and diffraction of 
energy may be neglected. Hence, since the flux of energy across any 
section of the ray tube is proportional to the square of the amplitude 
multiplied by the cross-sectional area of the tube, the amplitude varies 
with distance s along the tube like A~'*(s), where A(s) is proportional to 
the cross-sectional area at the point s. 

But, even when there is no reflection or diffraction of energy, the 
dissipation of energy by the shock wave (which in reality replaces the 
wavefront) and the related distortion of the wave profile behind the shock 
due to non-linear effects cannot be ignored. ‘Thus, even for weak shocks, 
the result that the shock strength varies with s like A~!*(s) requires 
modification, just as in the symmetric problems of the original theory, 
It should be stressed that this inaccuracy is a failure of the linear theory of 
sound and is not introduced by the approximations of geometrical acoustics. 
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Now, although its prediction of the shock strength is incorrect, 
geometrical acoustics provides the key to the solution of these more general 
shock problems: we assume for them also that the propagation of the distur- 
bance down each ray tube may be treated separately. ‘This gives a two variable 
problem depending on time ¢ and distance s, and it can be solved by precisely 
those methods which were developed in the original theory. The other 
variables in the problem appear only as parameters in the function A(s) 
and in the function which specifies the initial wave profile for each tube. 

In the improved theory, any point of a shock moves with the speed 
appropriate to the strength of the shock at that point. ‘Thus, even for 
propagation into a uniform medium, if the strength varies along the shock 
there will be a tendency for the shock to be refracted away from the wavefront 
positions given by the linear theory (which assumes the uniform speed of 
propagation a,). As a consequence, the true orthogonal trajectories of the 
shock positions will curve away from the straight rays of the linear theory. 
In principle, therefore, the ray tubes need modification at the same time 
as the modification to the law of propagation in each tube. However, 
unless the strength varies very rapidly along the shock, the effect of the 
curvature of the rays is relatively small. ‘he displacement of the ray 
from its linear position may become large as s-- x, but the displacement 
remains small compared to s, and the total angle turned by the ray 1s small. 
Since we are most interested in the directional distribution of shock strength 
for given s, this error may usually be neglected. 

Nevertheless, situations do arise in which the shock strength varies 
rapidly along the shock. An example of this occurs in the explosion problem 
when the surtace S is concave outwards in some region. For, then, the 
rays intersect, and since A~--0 at points of intersection, the linear theory 
predicts that the strength will become infinite. Of course, the linear theory 
breaks down completely and the consequences are shown in figure 1; 
AB is the initial surface, CD, EFG, HL/K are successive positions of the 
wavetront and the dotted line is the ‘caustic’, i.e. the envelope of the rays. 
Now, when the properties of a real shock are taken into account, this singular 
behaviour is avoided. In the concave part, due to the convergence of the 
flow, the shock increases in strength, and therefore moves faster than the 
neighbouring parts. ‘This effect smooths out the concavity and the rays 
curve away from each other, avoiding intersections in the neighbourhood 
of F. ‘The true type of behaviour is sketched in figure 2. In such a case 
the distortion of the rays is crucial. By making simplifying assumptions, 
this distortion can be studied to some extent mathematically and a theory 
can be given which covers the main features of figure 2. However, the 
theory forms part of a separate investigation which has more general 
applications; accordingly, it will be left to a later paper. In most of the 
problems treated here, this question does not arise, and it will be sufficient, 
for the overall picture, to treat only the propagation outside such singular 
regions, bearing in mind the qualitative description of the behaviour as 
represented in figure 2. Thus, the theory will be presented neglecting the 


deviation of the rays from their linear positions. 
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So far the basic notions have been described for the explosion problem, 
but the account applies in general to all the problems to be considered, 
provided that in supersonic flow problems the system of reference is so 
chosen that the air is at rest at infinity. In each case the linear theory 
based on the assumption of small disturbances) represents the shock as 


A C E i 
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Figure 1. 











Figure 2. 


a wavetront spreading out from the body or explosion with its amplitude 
varying as 4-'*(s), where A(s) is proportional to the area of the appropriate 
ray tube, and in deducing the correct results, the flow in each ray tube may 
be considered separately. Thus, when the general results for propagation 
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in a non-uniform tube have been established, the application to specific 
problems reduces to the determination of the area function A(s) and the 
initial wave profile for each ray tube. In those problems of supersonic 
How which reduce to steady flow problems when axes fixed with respect 
to the body are chosen, it is more convenient to treat the steady flow problem 
directly rather than to introduce the time as an additional variable. ‘The 
ideas developed above have similar interpretations in steady flow; these 
will be given and used in $$ 6 and 7 for examples from that field. 

Only propagation through a uniform medium will be considered, but 
it may be noted that this is not an essential restriction. ‘lhe theory of 
geometrical acoustics will always give the basic geometrical picture of 
wavefronts and rays, but when the sound speed is not constant the rays 
curve round, as the wavefront is refracted. The prediction given by the 
theory of sound for the amplitude of the distrubance may again be deduced 
from energy conservation along a narrow ray tube, but factors involving 
the density p, and sound speed a, in the undisturbed fluid are no longer 
constant and must be retained in the expression for the energy flux. For 
example, the amplitude of the pressure variation is proportional to 
(py) a@))'?A-'*. In principle, the corrected results for the propagation of the 


shocks can be deduced by the methods described below. 


In this section, a general account is given of the method for improving 
the linear theory of the propagation in individual ray tubes. ‘lhe method 
is developed from physical arguments; but, as a check, the results for 
the ultimate decay of the disturbance at large distances along the ray tube 
are also deduced mathematically. ‘The required non-linear features may 
be introduced by taking account of the progressive distortion of the wave 
profile due to the small variations in the values of the propagation speeds 
of the individual wavelets in the wave; each wavelet travels with the local 
speed of sound a relative to the fluid, and this is only approximately equal 
to a,. ‘The dissipation of energy at the shock is then incorporated auto- 
matically at a later stage, simply by applying the Rankine~Hugoniot shock 
conditions. 

The first step is to examine in more detail the results given by the theory 
of sound. For the most part, problems of shocks moving into undisturbed 
Huid are considered, and therefore it is appropriate to find the expressions 
for the flow quantities near the head of the wave. Now when the various 
problems are considered, it is found in general that in each ray tube the 
pressure increment p— /p, and particle velocity u are proportional to 


FU—s 4) 

\ A(s) 
where A(s) is the ray tube area and the function Ff’, which determines the 
detailed wave profile, depends on the initial conditions in the particular 
problem considered. ‘Thus, near the head of the wave, the amplitude ts 


(1) 
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correctly predicted by geometrical acoustics; however, the full solution 
has to be used in order to determine the function F. 
To include the non-linear distortion of the wave profile, (1) is modified 


to F(r) 
Feo) @) 

VA 
where 7(¢, s) is to be determined so that each wavelet specified by 7 = constant 
travels with the accurate speed a+u in place of a. Hence, 7 is to be 
determined from ds . 
(5) =at+u. (3) 

dt T ynstant 


Since the disturbance is assumed to be small, a—a, is proportional to 
p—p»; therefore, a+wu differs from a, by a multiple of F(z)/A. Then, 
to the same order of approximation, (3) may be written 


(5) 1 1 kF(r) (4) 
ds}, wo. 2 a A’ 
where & is a constant. ‘lhe arbitrary function of 7 which arises in the 


integration is fixed so that 7 takes the linear value t—s/a, near s = 0, in 
order that the wave profile agrees with (1) initially. ‘Then, (4) gives 
ds 


r= a, “RE (r) | JA rT. (: 


J 
— 


and the increasing divergence between the values of 7 and t—s/a, as s 
increases gives the progressive distortion of the wave profile. 

In all the cases considered, it is found that u = (p—pp)/po@, where 
py is the density of the undisturbed gas (the necessity for this result is seen 
below in equation (10)), and for a polytropic equation of state (p < p’”) 


a—aA y—1 p—Pp 





£8 ¢ 
Ay 2y i (0) 
Hence, if F(z) is chosen so that 
Po \ A’ 
then 
u 1 F(7) a—ady y—1 F(r) 
Se _ mat aapeeocaeas a ae *. 8 
a, y VA’ ay 2y A’ (8) 
and, in (4), 
~v/ ~~ 1 
c= aes - y 
2yay ) 


he improved solution to the problem is now given by (7) and (8), with 
z given implicitly as a function of tand s by (5). Inthe problems considered 
here, the wave is headed by a shock; therefore, to complete the solution, 
the position and strength of this shock must be determined in accordance 
with the Rankine-Hugoniot conditions. For a shock moving into un- 
disturbed gas, it may be shown that, to a first order of approximation, 


uy 1p, Po a, — 4 y¥—1 pi—DPo U iia 8.3 l p, —Po (10) 


: oa 9 
ay 4y Po 








dy Y Po a 2y Po 
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where U is the shock velocity and the subscript 1 denotes values behind 
the shock. ‘The first two of these equations are satisfied already in the 
solution (7) and (8) (even in the linear theory these conditions are satisfied 
by discontinuities in the flow quantities), and the third condition determines 
the shock. It may also be noted that the shock conditions show that the 
entropy jump at the shock is third order in the shock strength, which 
explains why the entropy changes can be neglected in the flow behind the 
shock; this point is discussed in detail in Whitham (1952). 

Now, the shock is intimately connected with the distortion of the wave 
protile. Wavelets are continually fed into the shock, where their energy is 
dissipated. Without the shock, the wave profile would * break’, exactly 
as in the well-known theory of plane waves of finit? amplitude, due to the 
higher propagation speed of the wavelets in the regéons of higher pressure. 
Thus, the compression regions of the initial wave profile shown in figure 3(1) 
would eventually break, as shown by the full line curve in figure 3(ii). “This 





) (11) 


Figure 3. 


leads to a solution which is physically unreal since it predicts more than 
one value of the pressure at certain points. However, shocks cut out the 
overlapping wavelets as shown by the dotted lines in figure 3(ii); for example, 
at the head of the wave the pressure jumps almost discontinuously from A to B. 
It is clear that the position and strength of the shock are determined if 
the value of + corresponding to the flow just behind the shock is known; 
that is, if the value of + at the point B on the wave profile is known. — For, 
denoting this value of + by 7T(s), the time ¢ at which the shock reaches the 
position s is found by substituting 7 = 7(s) in (5), and the pressure, particle 
velocity ete., are found by making this substitution in(7)and (8). However. 
these results must conform with the shock conditions (10), and this deter- 
mines the function 7(s). It follows from (5) that 

| dt 1 kF(T) ton a ee 
-—~—+-~ +(1-kF(T) | —]—, 

l ds Qy \ A ( Jor a ds 
and, according to (1), (again retaining only the first order term), this must 


be equal to 


(1 i ohh) _ t my) 


ji wa 
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Hence, 








vf {? ds\dT 1 kF(T) dT 
kF(T) | a + a ( me 
IgpVA/ds 2 WA ds 
\fter multiplication by 2 F(T), this equation integrates to 
i 
2| WT’) dT’ 
ae (11) 
| evs RAFAT) 6 OC’ 
issuming for definiteness that 7 = 0 at s = 0. 

It is found in general that the positive phase of the wave (p— py > Y) 
in the region immediately behind the shock, is followed by a negative phase 
(p—py <0). When this is the case, P(r) has a zero for some finite value 
7, of 7 and the results for the ultimate decay of the shock take a simple form. 
For, assuming that | 4 !'*ds-~- ~ as s-~ %, (11) shows that 7(s) - 7) as 


s-- oo; and, in fact, (11) may be approximated for large s as 





2,;7T var rs ds \-12 
At it fs ee nt T'\ yf om Be aie ¥ p, 
MT) ~ 53) PU) a 1 itz (12) 
Hence, from (7) and (9), the value of the pressure behind the shock is given 
by . T, 12 8 F 12 
; Pi-Po _- tao | F(T’) aT’ 1A | ds on (13) 
Po pee l 0 ov A 


and u, and (a, —4,)/a) are proportional to this. From (5), the position of 
the shock at time ¢ is given by 
s oe en ee 
t atin | Mey. <] 


ay Y% <0 o\ A} 








+ Ty. (14) 


‘Thus, at s the shock is ahead of the wavelet t = 7,+5/a,, on which the 


pressure is zero, by an amount 
I(s)o<| ——> . (15) 


‘lhe quantity / 1s a measure of the length of the wave and (15) shows how it 
ultimately increases with s (in contrast to the constant length in the linear 
theory). Ihe pressure distribution within the wave also takes a simple 
form. For, at all points 7 is near 7,; therefore, equation (5) which deter- 
mines 7(t,s) may be approximated as 








s . * ds Aa 
t= Pa —kF(r) i a i | 
Hence, using this value for F(z), (7) becomes 

th. 1 Ss (16) 

Po kV/A oVA 
‘Thus, at any point s, the pressure falls nearly with time and the rate of 
fall is d(p—p i, *s ds )-1 = 
(PoP) = - ay! va |’ . (17) 
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It is of interest to note that (17) depends only on the distance s and the 
properties ot the fluid; it is independent of the initial wave form given 
by F(z). 

An essential point must be made in connection with these results. 
‘The expression (1) was originally introduced as the appropriate approxima- 
tion to the solution given by the theory of sound for the flow quantities neat 
the head of the wave. In the improved theory, wavelets are continually 
fed into the shock; hence, values of the flow quantities at points away 
from the head of the wave in the linear theory must also be considered. 
However, the wavelets are fed into the shock at a relatively slow rate and by 
the time this effect has become appreciable the wave has travelled a relatively 
large distance. It is found that (1) again applies for large s even when 
t—S/d, is not small (in fact the precise condition for the approximation (1) 
is usually that (a@,t — s) sshould be small); hence its use is justified throughout 
the motion of the shock. 

The results given by (13) and (15) tor the important cases of plane, 
cvlindrical and spherical waves may be singled out for special note. For 


plane waves, 
Ad = constant, (P1—Po)/Po K 5**, poe, 18) 
tor evlindrical waves, 
Acs, (Pi—Po) Po % 34, loc sh 4, (19) 
and for spherical waves, 
As’, (P1 —Po)/Po & s(log s)-¥?, / x (log s)!?. (20) 


Using the simplified results (13) and (15), it is possible to see the 
significance of the existence of a zero of F(z). First, consider the outward 
flux of mass at the point s as the wave passes that point. Clearly, the 
Hux due to the positive phase of the wave, in which F(z) > 0, is proportional 
to po, Al; according to (13) and (15), this quantity varies as 1/A. But 
apart from the case of plane waves for which A is constant, it is generally 
true that 4 4d -- cass-- a. Inthe latter case, therefore, it is impossible 
for the wave to have only a positive phase, and so it must be followed by 
a negative phase; hence F(z) has a zero. Furthermore, in any realistic 
case, the fluid returns to a state of rest after the whole wave passes to that 
there is a recompression region after the negative phase, and a second 
shock must ultimately be formed (due to the distortion and ‘breaking’ 
ot the wave profile). Since the positive and negative phases of the wave 


must eventually produce equal and opposite contributions to the mass 
flux, the wave profile will tend to a symmetrical N-wave in which the two 
shocks have equal strengths and the rate of fall of p between the shocks is 
given by (17). The formation and subsequent motion of the second shock 
can be described in detail by the present theory, and the tendency to form 

symmetrical \-wave is confirmed. A full account of this is given in 
Whitham (1952) for the special cases of steady supersonic flow past an 
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axisymmetric body and unsteady plane waves; the extensions to the more 
general case considered here are trivial. 

At this stage, having described the main structure of the theory, it 
is suitable to reconsider the physical assumptions that have been made, 
and to give mathematical checks where possible. ‘There are two points 
which require further verification. ‘The first concerns tlie energy dissipation 
by the shock; this was not introduced explicitly and the question arises of 
whether it is correctly accounted for by making use of the shock conditions. 
Now, it is easily verified from the results in (13) and (15) that the wave 


loses energy at exactly the same rate as the shocks dissipate energy. ‘The 
full details will not be given, but we may note that the two quantities have 
the same dependence ons. ‘lhe energy carried by the wave ts proportional 
poaial( > Pe)’; 
Po 
from (13) and (15), this quantity varies with s like 
$ ds )-12 
} oVA 


he rate of dissipation of energy by the shock is proportional to 


pyas( PsP A, (22) 


0 


the increase of entropy at the shock being proportional to the cube of the 
shock strength. From (13) we see that (22) also varies with s like (21); 
when the constants of proportionality are included, exact agreement is 
found. 

The second point which may be considered further concerns the basic 
step which introduces the non-linear distortion of the wave profile, i.e. the 
introduction of a+u for the signal speed of individual wavelets. Of course, 
this rests on sound physical argument, and indeed is exactly what is found 
mathematically in the well-known theory of simple waves in one-dimensional 
unsteady flow. Nevertheless, further mathematical justification is desirable. 
To provide this the simplified form of the results for large s ((13), (15), (17) 
etc.) will be established by an alternative method, which proceeds directly 
from the non-linear equations of motion. In fact this method gives the 
simplest derivation of these results. 

Since the entropy jump at the shock is of third order in the strength, 
entropy changes may be ignored in the flow behind the shock. (This is 
borne out by the above considerations of energy balance.) ‘The equations 


U2 
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for one-dimensional flow in a tube of cross-sectional area A(s) may therefore 


be taken as 


u Cu 1 op (23) 
= rea: = a 
ct cs p cs 
Pp oa p". 
Introducing the sound speed a, equations for u and a may be written 
Ca ca =6y—l feu A’ 
— +y— + —— aj — + — uw] = 0, (24) 
ct cs 2 cs A 
Cu ou z Ca ia 
— +u— + —a- 0) (25) 
ct cs aan os 


With the N-wave in mind, let us now consider the possibility of the following 


solution in series: 


4 A) § \« 
u o(9)( he -) : v4(s)(t d ies 4 Ee. (26) 
dy 2 dy 
a= a,4 dy(s)(t Sie a, . by(ay(t = ~ ) ee (27) 
Ay : ay 


(The constant 7, is introduced in the measurement of ¢, in order to agree 
with the previous results.) Substituting these expressions in (24) and (25) 
and equating coefhcients of (t— T,, —s/ay), it is easily found that 


y-1 
b,(s) : ae w'y(S) (28) 
and 
de yt] ve = a v1 (29) 
ds 2 a 2A 


Equation (29) may be written 


=| l y+1 1 Q- 
aoe * a ca 


hence, 


—, (30) 


If only the first order terms are retained in (26) and (27), we have exactly 
the solution found earlier; for example, 


p Po 2y a= A 2y t re ZT, — $/Qo ( 4 1 ) 
— — ——— = — —— gq, ——_—_————— 3 
. y—-l a y+1° salt : 
P ‘ , ‘ \ A ds/4 A 


in agreement with (16) and (17). Moreover, if the head shock is specitied 


by 


i ihe oA, (32) 


ay 











On the propagation of weak shock waves 301 
so that its velocity is 
U ( ] dl\-) - dl 
=(—-— = Ay—-A >, 
a, ds a 


then the shock condition (10) shows (using (31)) that 


dl ik (" ds I 
ds 2/Al !9 VA 
‘The solution of this equation is 
re ds \12 
(s) " | | 7A , 


in agreement with (15). ‘he shock strength is obtained by substituting 
(32) in (31). 

In this derivation, questions of convergence have been ignored, but it 
is expected that the ratio of successive coefficients in (26) and (27) decrease 
essentially like 1/s and that the series are at least asymptotically valid when 
(t — 7, —s/a,) is small compared with s/a,); a(t — T, —s/ay)/s is less than 1/s, 
and in all the problems considered //s tends to zero as s—~ 0. 

‘The above arguments verify the results obtained for large s quite generally. 
Other checks on the theory will be made in specific problems by comparing 
the predictions of the initial shock strengths with those obtained by other 


methods. 


3. SPHERICAL WAVES 
Before considering the unsymmetrical problems which are the main 
subject of this paper, it seems worthwhile, in order to illustrate the general 
theory given in § 2, to include the simple example of spherical waves. 
First we note that each elementary ray tube is a cone with vertex at 
the centre of symmetry; therefore, if s is measured from the centre, we may 
choose 


A(s) = 82. (33) 


Secondly, we must consider the linear theory and verify that the results 
quoted in §2 as being typical of all problems, are found in this case. In 
the theory of sound, the flow quantities may be deduced from a velocity 
potential ¢ which satisfies the wave equation; for spherically symmetric 
waves the appropriate solution is 


‘os f(t—s a) 
Ss ; 
Now, u = 66/0s, p— py = —potd/ct; hence, 
Pp — Po 2» ¥ f(t -$/) 
Po a s : 


1 f’(t—s/a,) f(t — §/Ay) 


a \) $$” 
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In general the wavefront will be given by t—(s—s,)/a) = 0, where s = s, 
is its position at ¢=Q; it is convenient, therefore, to introduce 
7 =t—(s—s,)/a,. Then, if F(r) is defined by 

: "7a Eigen ; 

}(r) = | F(z’) dz’, 


the expressions for (p— py) Py and u become 


) F (7) 
P—Py | aor 
Pp S 
u 1 F(r "Sag (5 eee eae pe 
) — | F(t’) dr’ . (33) 
a, Y r ys” J 


When 4,7/s 1s small, the second term in (35) may be neglected in comparison 
with the first, and we see that the results quoted in (7) and (8) are borne 
out by this example. 

The form of the additional term in (35) is of interest in view of earlier 
remarks. For, in any physically realistic case, p returns to p, and uv returns 
to zero after the wave has passed. But we see that when F returns to zero 
after a positive phase, w still ditfers trom zero by the (smaller) second term 
in (35); in order to reduce both p to p, and u to zero there must be a further 


negative phase so that | F(z’) d>’ -0. This is an alternative argument 


for the existence of a negative phase to the wave. However, it is closely 
connected with the mass flow argument. For, to produce the large outward 
mass flux (proportional to y A = s) in the positive part of the wave, there 


must be a continual net forward transfer of mass across the wavelet t = 7, 
which separates the two phases. ‘This is represented in (34) and (35) by 
the fact that onz = 7), where p = /p,, wis zero only to the first approximation, 


and the second term in (35) gives the required flux. 

In any specific problem, F(z) is determined by the boundary or initial 
conditions. ‘These usually take the form of a prescribed value of u on some 
surface s = R(t) (for example, a particle path may be prescribed), and F(z) 


is obtained from (35) by solving the first order linear equation for | F(+’) dr’. 


It may be noted that in most cases the surface R(t) will not be in the region 
where the second term of (35) can be neglected and therefore the theory of 
geometrical acoustics cannot be used throughout; in fact, if R(t) is small 
the second term is dominant. 

In the improved theory, there is little to add to the results in $2; one 
point, however, requires care. Let us assume that the disturbance is 
generated at the surface of the sphere s = R(t). If the initial radius R(O) 
is not zero and R(t) is approximately equal to R(O), equations like (5) and 
(11) require only slight modification to take account of the fact that the waves 
start at s = R(Q) rather than ats = 0. For example, (5) becomes 


s— R(O ; ds 
da = ) kF(r) aay tae (36) 
ay RO) SN A 
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But, if R(O) = 0, more care is required since with 4 = s*, [ds/\/A is not 
convergent at s = 0. However, the lower limit is chosen so that 7 agrees 
approximately with its linear value for the initial propagation of the wave 
form; hence, we may, with greater accuracy, replace R(0) in (36) by R(7), 
where R(r) is the value of R when the wavelet labelled by 7 is generated at 


the sphere. ‘Then, using A(s) = s*, 7 is determined from 
s— R(r) s 
f= —— - kF(r)log —— +7. (37 
a, R(r) en 


Similarly, the relation for 7(s), the value of 7 at the shock, becomes 


1 
>| F(T’)dT 
: ‘ a ee 38 
R(T) k FUT)" ia 


§ ys 
log — 

‘To provide a check on the theory, we may consider the shock produced 
by a sphere expanding at a uniform rate, since this special case has been 
solved using a different method by Lighthill (1948). If the radius of the 
sphere at time ¢ is R(t) = ma,t, the boundary condition is 


u = May on R = mat. 


From (35), this determines F(r) for small values of m to be 

F (7) = 2ymbayr. 
Ihus, since F(7) is continuous at 7 = 0, there is no pressure jump according 
to linear theory. But, in the improved theory a shock is predicted, as 
required, although its strength is extremely small. Equation (38) gives 


; 7 l l 


~~ = ene 
ma, 1 k2ym?a, (y+ 1)m? 


Therefore, at the shock, 


ce of By & ; 
Pi~Po = F(T) 2ymre Ly-3 


, 


Po Ss 
and the shock velocity is given by 
U y+ ee 
ee able nS 
a 2 


0 
In these expressions, the factors multiplying the exponential are suspect, 
since the error terms in the exponential may well be more important. ‘The 
most that can be said with certainty is that 


l l 
log (= -1)~ a. 
Ay (y + 1)m* 


This is the result obtained by Lighthill. 


+. UNSYMMETRICAL EXPLOSIONS 


In this section the theory is applied in detail to the explosion model 
described in §1. We consider a high pressure region V of arbitrary shape 
in which the gas is at rest at a pressure p,+P, and which is surrounded 











304 G. B. Whitham 


by an infinite expanse of gas at rest with pressure p,; at time ¢ = 0, the high 
pressure region is released. Since only weak shocks are considered in this 
theory, it is assumed that P/p, is small. Furthermore it will be assumed 
that the region Vis convex; the difficulties which arise when this is not the 
case will be noted later. 

‘The rays in this case are normals to the surtace of |’, and the expression 
for A(s) is easily found in terms of the curvature of the surface at the foot 
of the ray. If R, and R, denote the principal radii of curvature at any point 
P ot the surface, the radii of curvature of the wavefront at a distance s out 
along the ray from P are R,+s and R,+s. ‘Therefore, by considering the 
area of a small curvilinear rectangle formed by the principal curves on the 
surface, it is seen that the area of any small element of the wavefront Is 
proportional to (R,+s)(R,+s). ‘Thus, we may take 

A(s) = (R, - s)(R, +S). (39) 

We next verify the form of the solution (7) and (8) and determine the 
function F trom the full linear theory. In the linearized formulation of 
the problem the velocity potential (x, Vv, 3, ¢) satisfies the wave equation. 


V°¢ = = 
and the initial conditions 


@ (), = = Po 
0) outside IV, 


where p, is the density of the gas. (‘The first condition arises since the 


velocity ¢ Vo is zero, and the second from the result that in the theory 
of sound the pressure excess p— p, is given by —ppjcd ct.) ‘This is a special 
case of a classical problem solved by Poisson and the derivation of the 
solution is given, for example, in Lamb (1932, Art. 287). “The solution for 


general initial values of d and ed ct Is 


Cd C 2 
d(x, y,2,t)=tM [3 ]- ~<tM, [9] ;, 
Tr Ge ct 


= 

where M, |] and M, [cd ct] represent the mean values of the initial 

values of d and cd/ct taken over the surface of a sphere with centre at 

(x, v, =) and radius af. ‘Thus, in the present case, 

Pr ge,9, 2,8) 
47p,a" t 


where g(x, v, 3, ¢) is the area of the part of the surface of the sphere which 


(40) 


lies inside J’. 

We now consider the variation of ¢ with ¢ at a point a distance s along 
the ray from the point P on the surface of )’. In particular, we consider 
the approximate form of the solution near the head of the wave and for large 
values of s. Clearly g remains zero until the radius at of the sphere reaches 
s, corresponding to the arrival of the wavefront att = s/a,; we setr = t—s a). 
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For small 7, the sphere will intersect V in a small curve surrounding P. 
If we choose coordinates (€, 7, ¢) with origin at P, ¢ measured along the 
ray at Pand and » in the principal directions at P, the equation of the surface 
of Vin the neighbourhood of P is approximately 


& 7° 
Cc=-7zy5>—-aye 4] 
aR, 2R» oO 
‘Lhe equation of the sphere is 
({—s +74 "° = (S+ay7)’, (42) 


and, assuming that a,7/s is small, its intersection with (41) satisties 


le ie l ef l : 
at ees eee ( ae ee (43) 
2a,T\ R, s Zaar\ Re § 
To the same approximation, the area of the surtace of the sphere bounded 
by this curve is equal to its projection on the plane ¢ = 0, i.e. it is the area of 
the ellipse (43). Therefore, 


R,R,s? 12 
9(x, y, 2, t) = 22agr, ——*+~_+——— . +4 
Bt. °0") (R, +5)(Ro +5) = 
Then, setting tf = s/a, in the denominator in (40), we have 
P RR; 12 : - 
ek aaa ncceaee 7, for small 7. (45) 


2p)|(R, +s)(Ro+5) 


It may be noted that p—p, = —p9cd ct jumps discontinuously from (0 to 
}P at s = 0, and this agrees with the well-known one-dimensional result. 
The solution is also required when 7 is not:small but s is large. For 
sufficiently large s, it is clear that the part of the sphere intercepted by V 
may be approximated by a plane perpendicular to the ray at P and at a 


distance a,7 along the ray inside 1’. ‘The area of this plane inside V is 
independent of s. ‘Therefore, again taking ¢ = s/a, in (40), we have 
P f(z) 
d= — d 46 
47pjQy Ss ia 


where f(z) is the cross-sectional area of V’ at a distance a,7 along the inside 
normal from the point P, the section being taken perpendicular to the normal. 
Only the dependence of f on 7 is shown, but it must be remembered that it 
also varies from ray to ray. ‘lo the same order of approximation, we may 
write 

P f(z) 47 
~ Tapa, VR + (RH) (#7) 
which brings the result into the same form as (45). In fact (45) may be 
included in (47), provided that f(r) ~ 27a 7\/(R,R.) tor small 7. But this 
is so; the intersection of ¢ = —a,7 with (41) is an ellipse with semiaxes 
\/(2a97R,), 1/(2a,7R,), and the result follows. Hence, (47) applies both 
for small 7 and large s, and the appropriate requirement is that @y7 s is small. 
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As noted in (38), A(s) <(R,+s)(R, +5); thus, (47) is of the form 
predicted by geometrical acoustics. We have 


6) Ch F : 
p- Po ee Fv oo. (48) 
. Po Po ct \ 4 
wnere > 
str) 
F(r) = —, (49) 
Dy 47a, 


in accordance with (7). Moreover, the particle velocity, o¢/és, along the 
ray agrees with (8). Thus, all the results quoted in the general theory 
of §$1 & 2 are illustrated in this example. ‘The non-linear theory derived 
in $2 applies directly to this problem, with the functions A(s) and F(r) 
for each ray given by (38) and (49), respectively. Inthe theory, | 4-!?ds 


is required, and we may note that it is given by 


= Ying 4 » (Ri +s) + V(Re+5) (50) 
\VA ia \ R,+ VR, : ; 
The shock strength is given as a function of s by (48) with 7 = 7(s) as 
determined by (11). 
Since F(7) tends to the finite value | Pp, as 7 -- 0 the initial strength 


of the shock is }P/p, as in the corresponding problem of one-dimensional 
flow; the compression wave increasing p from p, to py+}P travels out 
into the undisturbed medium and an expansion wave reducing the pressure 
from p,+ P to p,+5P travels into the high pressure region. ‘Then, for 
large values of s, the simple law (13) applies, and 4 is proportional to s* 
so that the variation of shock strength with s is just as for a spherical wave. 
But, the constant of proportionality in (13) varies with the direction of the 
ray. Itisseen from (49) that the zero 7, of F(z) corresponds to the maximum 
cross-sectional area of V perpendicular to the ray considered. The function 
f(z) increases from zero to a maximum at 7 = TJ, and then decreases to zero; 
hence F(r) takes both positive and negative values and again 


} | F (7) dz 0). 
From (49), 
T P 
| F(r)dr z 
0 7) tna,p,. 
hence, (13) gives 
ne 4 ve 
—— ~ - F0T s) s-, log ——__——_—_——; > i (51) 
Pi y+1 7 pp (y Ri +4 R,)* 


for large s. ‘he shock strength at a tixed distance is proportional to the 
square root of the excess pressure and to the square root of the maximum 


cross-sectional area of V for that direction. ‘Thus the shock is strongest 
in the directions for which the projected area of V is greatest. Even for 
strong explosions, we may expect the predictions of the theory to be 
qualitatively correct, and the directional variation of shock strength to be 
correlated with the projected areas of V. 
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It is perhaps worth noting that for the special case in which V is a sphere 
of radius R,, (47) is the exact solution for all sand 7.. In this case 
f(r) = n(2R, — agr)ay7, 0< a7 < ZK, 
‘so that 


> 
2p, 
‘Thus, the disturbance is an \-wave from its inception, not only at distances 
sufficiently large for (16) to apply. 

If the boundary of V includes a region which is concave outwards, 
the rays will be distorted as shown in figure 2; but, after this region has 
been left behind, the rays will diverge, and ultimately the shock will follow 
the usual s (log s)'? law of decay. Equation (51) will probably still give 
a reasonably accurate form for the constant of proportionality. 

It should be noted that even if the boundary is not concave but is plane 


Fir) Ry — a7), O<a,r <= ZR,. 


over some small region distortion of the rays is still required to obtain an 
accurate result. For the ray tube area A(s) is constant on the rays from 
a plane region, and the propagation along these rays is initially similar to 
the case of plane waves. Hence the shock will not decay with s as rapidly 
as on the neighbouring diverging rays. Now this is the essential point: 
the dependence on s follows a ditferent law so that if unchecked the relative 
ditference in shock strength becomes large as s - «. But, starting first 
at the edge of the plane region, the effect of the large difference in shock 
strength will be to curve the rays outwards, and, eventually, they all diverge 
giving decay of the spherical type. It is interesting to observe that this 
view is confirmed by the expressions found for ¢: (45) applies for small 7 
and gives the plane wave formula as R, and R, tend to infinity, but (46) 
still holds for large s and gives the spherical form for ¢. 

When concave regions of the boundary V are admitted, the function 
f(z) may have more than one maximum for some directions. If this is so, 
there will be compressian regions in which F’(7) — 0 in addition to the two 
main compressions, and this leads to the formation of additional shocks. 
These may be determined in the same way as the multiple tail shocks in 
the problem of the supersonic projectile (Whitham 1952), but ultimately 


they run into one of the two main shocks. 


5. SUPERSONIC BANGS 


\n important application of the general theory is to the determination 
of the shocks produced by a body in non-uniform supersonic flight. ‘This 
application has been developed in detail by P. Sambasiva Rao (1956a & b), 
and the main results are quoted here for completeness. 

Near the body, just as in the case of uniform supersonic motion, the 
wavefront forms a cone about the direction of flight with semi-angle equal 
to the Mach angle sin-! 1/4, where M is the Mach number of the body 
at that point. Hence, the rays from any point of the flight path initially 
make an angle cos '1 .W with the direction of motion. But if the medium 
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is assumed to be uniform, the rays must remain straight. Hence, we have 
the typical pattern shown in figure 4 for an accelerating source. If the 
velocity were constant, the area of each ray tube would be increasing 
proportional to s, due to the cylindrical spreading of the wavefront away from 





Figure 4+ 

the axis. But if the body accelerates, cos"! 1/M increases and the rays in 
any meridian plane converge; conversely, if the body decelerates, the rays 
diverge. It is easily shown that the modification due to this extra 
convergence or divergence is included by taking 


5 Pe aes =9 

- I(s) (1 i} (52) 
where a,(M?—1) 

VW 'dM/dt’ 


Vf and dM dt referring to the Mach number of the body when it is at th¢ 


(33) 


foot of the ray concerned. For acceleration, A > 0, and we see that 4 
vanishes when s = A; at this point the wavefront is cusped. In figure 4, 
the rays shown as full lines ‘carry’ the forward part of the wavefront and 
s - Aon each of these rays; the rays shown as broken lines ‘ carry’ the rear 
part of the wavetront and on each of these s has exceeded A. For a curved 
Hight path Rao (1956) shows that 4(s) 1s still given by (52) but A is modified. 
he denominator of A is the component along the ray of the acceleration 
of the body; for a curved path this includes a term due to the transverse 
acceleration of the body. 

If the acceleration of the body is relatively small (the change of velocity 
in flying a body length is small, say) we might cxpect the initial wave profile 
to be the same as in the case of uniform motion; Rao’s detailed investigations 
contirm this. ‘lherefore, the /-function is the same as in Whitham (1952). 
Introducing a constant of proportionality in order that F conforms to the 


definition (7), 1t is found that 


’ l Ve 12 | Ur S"(€) dé 4 
F(t) = —\ ee er et (54) 
AMP-1){ te}, WUr—é 
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where S(€) is the cross-sectional area of the body at a distance € from the 
nose, and as usual, 7 measures the time after the wavefront passed the point. 

With A(s) given by (52) and F(r) by (54), the general results of §2 may 
be used. For s <A, A(s) ~ s and the shock spreads away from the axis 
like a cylindrical shock. If the velocity of the body is uniform, A is infinite 
and this behaviour extends right to infinity with the ultimate decay in shock 
strength like s-34. But if A is finite, we see that for sufficiently large s, 
A «x s* so that the decay is eventually like a spherical shock with strength 
proportional to s-'(logs)'?. It should be remembered, however, that 
when A ~ 0, a cusp will have intervened before this region is reached, and 
in the neighbourhood of the cusp the simple theory of geometrical acoustics 
breaks down. Hence, we are assuming that beyond the cusp the theory 
may again be applied, at least so far as the main features are concerned. 
The question of what happens near a cusp and what effects it has on the 
results is discussed in more detail in Rao (1956 a & b). 

Here, only the highlights of the theory of supersonic bangs have been 
mentioned in order to show the generality of the theory developed in §2; 
for the details and practical predictions, reference should be made to Rao’s 


papers. 


6. STEADY SUPERSONIC FLOW: AN EXAMPLE IN CONE THEORY 

‘The theory developed in the first part of the paper can be interpreted 
and applied in certain problems of steady supersonic flow past unsymmetrical 
bodies. If we introduce cylindrical polar co-ordinates (7, @,) and let L 
be the velocity of the undisturbed stream, then the steady flow problem is 
analogous to an unsteady flow in the (r,@) plane with x U playing the role 
of time. ‘Thus, for a pointed body, the analogy in unsteady flow would be 
to the disturbance produced by a solid cylinder which starts with radius 
zero and expands with arbitrary shape. If the body is swept behind the 
Mach cone through the nose, the corresponding rate of expansion of the 
cylinder is subsonic. In that case, the rays would be straight lines through 
the origin r = 0, and the linear theory would predict amplitudes proportional 
tor'?. Hence, for the steady problem, the analogue to treating the propaga- 
tion in each ray tube separately is to consider the flow in each meridian 
plane 6 = constant separately. Near the Mach cone, which is the analogue 
of the wavefront, the amplitude of the disturbance will be proportional to 
r- 12, ‘Thus any dependence on @ must be introduced by the profile function F. 

If we consider a body which is everywhere slender (i.e. its siope in the 
stream direction is always small) and whose cross-section has no regions 
of abnormally large curvature, the linear theory developed by Ward (1949) 
shows that, although the flow near the body varies with 6, the F-function 
yiving the flow near the Mach cone is independent of 6. Hence, the shock 
is the same as for a body of revolution with the same distribution of cross- 
sectional area. In order to obtain an example involving an unsymmetrical 
shock, we must relax these conditions on the body shape. 
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We consider the flow past a flat plate delta wing at incidence, the edges 
of the wing being swept back well behind the Mach cone. ‘This is a problem 
in which there is no fundamental length and hence the velocity components 
must be equal to the main stream velocity U multiplied by functions of 
r/x and @; the velocity potential has an additional factor » since it has an 
extra dimension of length. ‘This is a so-called cone field problem, and the 
linearized solution is known. Furthermore, Lighthill (1949) has shown 
how the shock strengths may also be found in such cone field problems. 
Here, the results are derived by the general theory which does not rely on 
the special properties of cone fields; the agreement with Lighthill’s result 
gives independent confirmation of the assumptions of the theory. 

A full account of the linear theory is given by Goldstein and Ward 
(1950). ‘The velocity potential is taken as U(x+¢) where 4 = xf(Br/x, 9), 
B = ,/(M?—1), and it may be shown for the present problem that, neat 
the Mach cone x = Br, 


j~ 2 9(8)( 1 — =). (55) 


The function g(#) introduced here ts identical with the 4(#) used by Lighthill, 
and is given by 
Lay 2sind P 

. = (56) 


(9) ni ne 
BI t, cos 6)?.*(1 t, cos #)3* 


where @ is measured from the plane of the wing, « is the angle of incidence 
of the wing, cot! t,/B and cot-!t,/B are the angles made with the ¥ axis 
by the leading edges of the wing, and 
es (fo ty) = — 
2(1 + t))!2(1 —t,)¥7{2E(k) —(1 — Rk?) K(R)}? 
(1-4 ty )(1 —t,) 
(1-4)(1+%)’ 
where E(k) and A(R) are the complete elliptic integrals. 
From (55), we see that 





L 


ke 


2 (8) Br)32 a 
Q=—;5 <a (* — ery"'* (I/ 
3 (Br)!? ) 
tor the flow near the Mach cone, and it may be noted that, as in previous 
examples, the crucial condition is that z/r should be small, where the 
characteristic variable tis now x— Br. Moreover, we see that the azimuthal 
component of the perturbation velocity, Ud,/r, is of order 7/r times the other 
components U¢,U¢,. ‘Thus, toa first approximation, the flow is in meridian 
planes and @ plays the role of a parameter. ‘This is in accordance with the 
general theory. From (57), we have 

je 

~ (Br)'* 
B'*9(0)(x« — Br)? 


yr} 2 
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or, introducing a standard notation, 


F(r) 
?, ae (2Br)!2” 
’ f (58) 
ae BF(r) 
ea (2Br)!2” 
where 
F(r) = g(0)2-¥ 722, (59) 


The improved non-linear theory is obtained by finding a more accurate 
relation for the characteristic variable 7(x,7r), in place of s=x—Br. This 
procedure is quite straightforward, and is the same as the corresponding 
step in the axisymmetrical problem (Whitham 1952) since @ appears only 
as a parameter. Briefly, it is as follows. ‘The characteristic direction 
at any point makes the local Mach angle » with the stream direction y. 
Hence, 7 is constant on each curve satisfying 


dx 


q, = Coty + Xp. 


‘The sound speed is given in terms of ¢ by Bernoulli’s equation; hence 
cot(u+ x) can be expressed in terms of ¢. Using ¢, = — B¢,, it is found 
that 7 is to be determined by 


dx : ; Y oy 4 
( ar iy ome ot ae 


to the first order in ¢. ‘Therefore, from (58), the relation for + is 


x = Br—k,F(r)r!? +7, (60) 
where 
(y+1)M! 
k, = “STERae (61) 


When F(z) < 0, the characteristics are diverging in an expansion wave ; 
but when F(r) > 0, they converge and form a shock. In the latter case the 
shock can be determined by giving the value 7(r) of 7 at points just behind 
the shock, and if the flow ahead of the shock is undisturbed, the condition 
1S 

1 
2| F(T")dT’ 
Pa. 62 
k,FAT) “— 
The derivation of this result from (60) is identical with that given in Whitham 
(1952); it also follows closely the derivation of (11) from (5). ‘The equation 
of the shock is then given by substituting st = 7(r) in (6). 

From Bernoulli’s equation, p — py = —p yU*d,, where py, py are the pressure 
and density in the main stream; hence, 

Ae y M24, sinnll d J (63) 
Po (2Br)! 


The shock strength is found by substituting + = 7(r) in this expression. 
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For the problem of the delta wing, F(r) <0 above the wing wher 
0 ..6< 7, and there is an expansion wave as expected. But below the 
wing F(r) >- 0 and there isashock. In this case, (62) gives 


so that 
Fi T) 3k, 27(6)) l 2 


Hence, the equation of the shock 1s 


3 (y+1)?M 
\ Bi T, - a g(A)r, (64) 
and its strength is 
) " M6 = 
= : iv(y + 1) 87(A). (63) 
0 ) 


‘hese expressions agree exactly with Lighthill’s results. 


‘THIN WINGS OF FINITE SPAN 
In the previous example only the profile function F depended on th 
orientation of the plane in which the flow was considered. We now turn 
to a problem in which the amplitude function (which was merely r~!? in 


the last example) also varies with orientation. This is the case in supersonic 


How past a wing with planform of the general shape shown in figure 35. 








Figure 5. 
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The linearized theory has been developed in detail (see, for example, the 
account in Ward (1955)) so that all the information required for applying 
the theory of this paper is obtainable. 

First we consider the geometry of the wavefront in order to deduce the 
set of planes in which the flow should be considered and to obtain the 
amplitude of the disturbances near the wavefront. We then turn to the 
linearized solution for the profile function F and, incidentally, for verification 
of the geometrical prediction of the amplitude. 

For simplicity, we consider oniy the problem in which the flow is 
symmetric relative to the mean plane of the wing; this is sometimes called 
the ‘thickness problem’. It is assumed that the boundary of the wing lies 
in the (x,y) plane with the v axis in the direction of the main stream, that 
the upper and lower surfaces of the wing are given by 

s= + Z(x, ¥), (66) 
and that the boundary of the wing planform is given by 
« = I(y). (67) 

The wavefront is the envelope of the Mach cones with vertices on the 
supersonic leading edge 4B of the wing, 4 and B being the points where the 
boundary makes the Mach angle p» with the main stream. If (/(y), 7) 1s 
a point on AB, the Mach cone is 

(x—I(n))? = B?(y— 7)? + B’s?. (68) 

At points on the envelope the derivative of (68) with respect to » is also 
satisfied, so that 

(~—/(n))l'() = B*(y—7). (69) 

The latter is the equation of a plane through (/(j), 7) parallel to the = axis 

and making an angle tan~!/'(7)/B? with the x axis. ‘This plane is shown 

is PNM in figure 6, where P is (I(n), 7) and PM is the intersection of the 


plane with the cone (68); PL isthe stream direction and PN isthe intersection 


of the plane (69) with s = 0. The wavefront is the surface generated by 
PM as P moves along AB. Now in applying the theory, the flow is con- 
sidered separately in each normal plane PLM. The coordinates 


corresponding to the ray coordinates of the unsteady flow problems are 


F.M. X 
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x,y and r, where 7 specifies the plane PLM and r measures distance from 
PL in this plane. Since LPN = tan-l'(n)/B and MPL = p = cotB, 
it follows that 6 = MLN is given by 





cos 6 = ad (70) 
Hence, the Cartesian coordinates (x, y, z) are related to (x, 7,1) by 
X= x, 
y=n+rcos8, (71) 
z=rsind, | 


The function corresponding to the area A(s) of a ray tube, is now the distance 
between two neighbouring normal planes PLM. If we consider the two 
planes corresponding to 7 and 7+6y and denote this distance by D(r)dn, 
the initial value of D is sin@; therefore, since the planes diverge with 
angle 40, D(r) is given by 


, dt 
D(r) = siné-r—-, 
dy 
\/ (B? - I'?) l" 





+ ar. 72 
B \/ (B? - 17) “a 
The amplitudes of the disturbances vary like D~-??. 

In order to find the solution near the wavefront in detail, we must turn 
to the full linear theory. But, we can already predict that the velocity 
potential ¢ will take the form 


. ~f7) - 
= TaBr+ (BT) BI’))’ (73) 





where 7, which measures distance behind the wavefront in a streamwise 
direction, is given by 


7 = x—I(y)—Br. (74) 


(Although only the dependence of f on 7 is shown explicitly it should be 
remembered that f is also a function of 7.) It is convenient to take the 
amplitude as given in (73) rather than as D~!?(r) itself, but it should be 
noted that this assumes /"(7)) + O at any point of the leading edge of the wing. 
The situation when /"(7) = 0 is very similar to the case of zero curvature 
of the surface in the problem of $4. In this case, the normal planes are 
parallel and the shock initially behaves as in the essentially two-dimensional 
problem of a swept back wing of constant section. However, as r increases, 
the normal planes will in reality curve away from each other to give the 
eventual decay typical of a finite body. 
For the symmetrical problem, the complete linearized solution is 


Z,(x', vy’) dx'dy' ” 
[(x— x’? — By —y’)? — B?s?] 





L ff u 
d(x, ¥, 2) = — . | | ; : (75) 
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where Z,(x, y) = 0Z(x, y)/ox, and the integration is over the part of the wing 
in the region x—x' > BY[(y-y'P +24]; (76) 
this region of integration is the intersection with the plane z = 0 of the 
upstream Mach cone from (x, y, 2). 

The approximate form of (75) is required both for small 7 and large r. 
These results turn out to be closely similar to the corresponding ones in 
$4, and the details of the check on (73) for small 7 is omitted. It is found 
that 

Ring A ae 

vl"(n) 

where e(7) = Z,(/(7), 7) is the slope of surface in the x-direction at the leading 

edge. ‘The approximation of (75) for large r is now obtained, and it gives 

a function f(r) which confirms (77). We set x’ = /(n)+z2, and y’ = 7+8, 

in the double integral (75); then, since the variations of x and § are small 
compared to 7, we have approximately 

(x — x’)? — By —y’)? — B22? 


9 2 l'r Pe 2,2 in 
= (Br+r+a)?-B (3 -B) — Br (1 ~ = 


2Br(r7 —x+1'(n)8). 


for small 7, (77) 





Thus, 








1 1;f 2, dads 
pad \ 2Br 7 ia V[r—(«—I'B))’ 
and the region of integration is bounded by the straight line «—/'B = 7 
which is parallel to the tangent to the leading edge at P and is at a distance 7 
downstream from P (see OR in figure 5). Hence, the integral in (89) is 
independent of r, and comparing (78) with (73) it is seen that 
1 ¢; 2, dzdB 
7) (r-@-TA) 

The expression (79) can be written in a form which agrees exactly with 
the corresponding result in flow past a body of revolution. ‘The region 
of integration is divided into elementary sections parallel to OR, then 
x—I'B is constant and equal to f, say, on any section. ‘The thickness of an 
elementary section is sinydt where y is the angle between OR and the 


(78) 





f(r) = 


} 


(79) 


stream direction; hence, for this slice, 
ff 2Zdad8 = S*(t)siny’ dt = S(t) dt 

where S*(t) is the cross-sectional area of the slice and S(t) denotes the 
projection of the area S*(t) perpendicular to the stream. ‘Thus, (79) 
becomes 

I)= 3) (80) 
and 

7) 
d= — Sr) (81) 


X 2 
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for large r. ‘This result is exactly the same as the velocity potential at 
distance r from the axisin flow past a body of revolution whose cross-sectional 
area at distance ¢ from the nose is S(t). ‘Thus, im each plane normal to the 
wavefront, the wing can be replaced by an equivalent body of revolution 
(as far as the flow at large distances is concerned). Of course, the equivalent 
body of revolution is not the same for each plane. 

The determination of the shock etc. is formally the same as in $6 with 
the slight modification that, in view of (73), r must be replaced by r+ r¢(7) 


where 
r() Br 


The velocity components and the pressure are given by 


- : 
F(r) ‘ BF(r) > 
fe" — EBs) 6°" VaBetry 
pp) _ _yMPF(x) 
Po \ [2B(r a ro)| ‘ 
where 
Ll. ¢t She 
F(z T ~ —_—_ 83 
(r) = f(r) te lea/r—D (53) 
The corrected expression to determine 7 1s 
y+1)M! 
\ 1 ”) t B — a F(r)\(r T ry 24 Ts (S4) 
At the shock, 7 = 7(r) where 
. (2py2 |. F(r) dt / 
(} r,)'2 Fate . (55) 


(y+1)M*” FT) 

The initial strength is given by setting 7’ = 0 in (82); hence, since 
F(0) f’(0) \ 2e(1) \ l” (see (77)), it 1s 
e(n)y.M? 


ie |i ba (86) 
\ | B?—-T%X)| : 


This is the same as the linear result, of course, and agrees also with the 
result for flow past a wing of uniform section and initial slope ¢, swept 
back through an angle tan“! /'(7). 

For large r, the law of decay is found by taking 7 near to the zero 7), of 


F(z), so that (85) gives the approximation 


ee (2B)3/2 -7 ; 1/2 ws 
ee doe | Br) dr} or, 87 
F(T) (+1 |, F( )d } (5/) 
Then, from (82), we obtain 
y  (2B)Ual (Te 12 
se (2B)™ | F(r)dr> . (88) 


Po (y4 1)!2 pt A 
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‘The shock decays in each normal plane exactly as in the meridian plane for 
a body of revolution having the appropriate F-curve, and the directional 
variation of strength arises only from the dependence of F on ». 
8. NOTE ON THE WAVE DRAG OF A FINITE WING 
The wave drag on a body can be calculated from the rate of dissipation 
of energy by the shocks. For a body of revolution the drag is then obtained 
in the form 
7p, U? | F(r) dr, (89) 
~ 0 
(see Whitham 1952), where F(r) is the F-function for the body. The drag 
on the wing can be calculated in a similar way. 
In the first instance, the contribution to the drag from the shock attached 
to the leading edge of the wing is 
21 | poTysD(r) drdn, (90) 
where p, and 7, are the density and temperature of the main stream, s is 
the entropy jump at the shock, D(r) is the function defined by (72) and 
1 <1 < Mp defines the leading edge of the wing. Now 


ee (y+1)U? Pi-Po\* 
6 “Damp )’ 


D(r) = —[r+ro(n)] - 
Hence substituting the value of (p,— fp )/p) from (82), the expression (90) 
becomes 





and 


~ A (y+1)M?! 





2 1 Sis) ks A Gy 4) 8 
te (2B) STN + ro ard, 
When r is eliminated by means of (85), this gives 
1 
pee (ee F(z) dz : 
py l . | 5 F%( I’) s ‘ |, \ y dTd# 
/ 0 0 dl ——Aa 
| F(T) 


“7 °T > 
ov? | (| F(T) a7) d§. (91) 


0 
Although the shocks from the trailing edge have not been treated here, 
n°?) 
it is reasonable to assume that they will contribute the terms | F°(T)dT7, 
. Ty 
since this is the contribution of the tail shock in the body of revolution case. 
If this is assumed the drag on the wing becomes 


pol? ( | FX) iT) dé. (92) 


/0 


~0U 
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On comparing this result with (89), it is seen that the drag on the wing is 
the mean of the drags on the equivalent bodies of revolution introduced in 
Section 7. Of course, the expression (92) for the drag must be equal to 
the value obtained by the more direct evaluation from the pressures on the 
wing surface, but for some purposes this may be a more useful form. It is 
certainly a much more significant form for the present theory. 
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Forces, moments, and added masses for Rankine bodies 
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SUMMARY 

The dynamical theory of the motion of a body through an 
inviscid and incompressible fluid has yielded three relations: a first, 
due to Kirchhoff, which expresses the force and moment acting on 
the body in terms of added masses; a second, initiated by Taylor, 
which expresses added masses in terms of singularities within the 
body; and a third, initiated by Lagally, which expresses the forces 
and moments in terms of these singularities. ‘The present investi- 
gation is concerned with generalizations of the ‘Taylor and Lagally 
theorems to include unsteady flow and arbitrary translational and 
rotational motion of the body, to present new and simple derivations 
of these theorems, and to compare the Kirchhoff and Lagally 
methods for obtaining forces and moments. In contrast with 
previous generalizations, the ‘Taylor theorem is derived when 
other boundaries are present; for the added-mass coefficients due 
to rotation alone, for which no relations were known, it is shown 
that these relations do not exist in general, although approximate 
ones are found for elongated bodies. ‘The derivation of the 
Lagally theorem leads to new terms, compact expressions for the 
force and moment, and the complete expression of the forces and 
moments in terms of singularities for elongated bodies. 


1. INTRODUCTION 

In the decade from 1920 to 1930 there were published by Lagally, Munk, 
and Taylor a number of hydrodynamic theorems concerning the added 
masses of bodies moving through an inviscid fluid and the forces and moments 
acting uponthem. These theorems enable the forces, moments, and added 
masses to be determined when the singularity distributions of sources, sinks, 
and doublets within the body, which may be considered to generate the 
potential flow about it, are known. Since, for the important class of elon- 
gated bodies, simple approximations to the singularity distributions are 
given directly in terms of the body shape, these theorems have furnished a 
powerful means of investigating the forces and moments acting on such 
bodies, especially near a free surface. Until recently, the scope of the afore- 
mentioned theorems has been limited: that of Taylor (1928) to only one 
kind of added-mass coefficient, and that of Lagally (1922) to steady flow only. 
Birkhoff (1953, p. 161) and Landweber (1956) have succeeded in generalizing 
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Taylor’s theorem to apply to all the added-mass coefficients except those for 
pure rotation, and Cummins (1953) in generalizing Lagally’s theorem to 
apply to cases of unsteady flows. 

The ‘Taylor theorem and its generalizations have heretofore been 
concerned with the added-mass coefficients corresponding to the motion of 
a single body in an otherwise undisturbed and unbounded fluid. ‘The 
extension of the theorem to the important cases where external singularities 
and other boundaries are present is highly desirable. 

Cummins was able to express the force and moment on a body in terms 
of the strengths of the singularities, except for one term in the expression 
for the moment—an integral over the surface of the body with integrand 
linear in the potential. As will be seen, this single unresolved term is 
intimately related to the missing relations in the generalization of ‘Taylor’s 
theorem. ‘The discovery of the latter would complete the generalizations 
of both the ‘Taylor and Lagally theorems. 

The present work, then, has several purposes. 

(1) The first is to extend the Taylor theorem to include cases with 
external singularities and boundaries and new results concerning 
the missing relations between added masses and singularities. ‘The 
latter will be derived for ellipsoids and for elongated bodies, but it 
will be proved that such relations do not exist in general. Also the 
opportunity will be taken to present, new, short, and simple proofs 
of the theorem for both two and three dimensional flows. 

(2) ‘he second is to consolidate and extend Cummins’ results and to 
present a simpler derivation of them. An important secondary 
motive for this part of the work is to popularize this powerful theorem, 
which in its present form has been applied only by Cummins himself. 

(3) The third is to examine the interconnections, if any, between the 
‘T'aylor theorem (relating added masses to singularities), the Lagally 
theorem (relating singularities to forces), and Kirchhoff’s equations 
of motion (relating forces to added masses). 


2. FORMULATION OF THE PROBLEM 
We are concerned with the interactions of a fluid with a rigid body 
moving through it. ‘he fluid is assumed to be incompressible and inviscid, 
the flow to be irrotational and, in general, unsteady. We shall suppose that 
the unsteadiness may be due to the time dependence of the linear or angular 
velocities of the body or to the presence of external boundaries or flow- 
producing mechanisms which may themselves be moving in an arbitrary 


manner. 

The flow may be considered to be generated by singularities, and the 
singularities considered will be isolated sources or sinks, doublets, and 
continuously distributed sources or sinks. Continuous distributions of 
doublets are excluded from consideration because they can be replaced by 
corresponding ones of sources and sinks. ‘The symbol m will denote the 
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strength of a source or sink, (a vector) will describe both the strength and 
the orientation of a doublet, and o will denote the volume or area density 
of continuously distributed sources or sinks. 

Cartesian coordinates (fixed in the body) x; (¢ = 1, 2,3) will be used, so 
that the position vector 7 is (x,, %2,¥%3) with magnitude r. The components 
of the velocity vector u of the origin of coordinates will be denoted by(u,, us, 
u;) and the components of the angular velocity w will be designated alter- 
natively by (w,, w,w ) or by (u4,u;,u,). ‘The surface of the body will be 
denoted by S, and those surrounding the singularities inside S by S’ collec- 
tively. The distance m normal to either S or S’ is always directed out of 
that portion of the fluid with which one is concerned. The direction of n 
will be indicated by the unit vector 7 with components , (i = 1, 2, 3), which 
are given by n, = ex,/dn. 





The kinematic boundary condition at a point on S is given by 
0p —- 
~~ =(ut+wxr).n=Uu,N,, (1) 

On 
in which ¢ is the velocity potential satisfying the Laplace equation, and the 
components of r x n are denoted by (m,,”;,,). Here the summation con- 
vention has been adapted and the Greek subscripts range from 1 to 6. For 
convenience of presentation, Greek subscripts will range from 1 to 6, whereas 
Latin ones will range over 1, 2, and 3 only, unless otherwise stated or when 
summation signs are expressly used. 

The velocity potential 6 may be considered to be composed of a part due 
to the motion of the body alone, when all other boundaries and external flow 
producing mechanisms are at rest, expressible in the form u, 4, after 
Kirchhoff, and a part ¢, due to the motions of the latter when the body is at 
rest, i. e. 


d= uy py + Po. (2) 

Then, on S, from (1) and (2), 
Od, _ Ody ae 2 
i ai) 2 _— oy = Ny. (3) 
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The boundary condition (1) may also be expressed in the alternative form 
U, My = (Uy + E54 Oj XN, 
where v, = —0¢/0x, is the -component of the velocity at a point of the fluid, 
and e,,, is zero if any two of the indices /, j, k are alike, and +1 or —1 ac- 
cording as the indices are in cyclic or anticyclic order. Hence, the boundary 
condition becomes 


Vm, = 0, Vi = Uj, +Uj, +i, WX, - (4) 


Since the coordinate axes are in motion, the Bernoulli equation for the 
pressure is, from Lamb (1932, p. 20), 
p Om : . ; "= 
- —W, W = 3(u,uj,+v,v,;)—v,(U, + €),) Wy X), (5) 
where p is the density of the fluid. 
It will also be convenient to have the expression for the kinetic energy 
T, of the displaced fluid, considered as a rigid body. Expansion of the 
B s ; 
integrand in 
2T p p | (u TE; 5p W; X;.)(U, T Ej, W) X)a7, 
where dr denotes an element of volume of the body, readily yields the 
quadratic form 


Uh 
2T, Bg uy ug, (6) 
B 8 Ba, : B,, = B6,; 9 B; j — Be jj.X ’ } 
, (7) 
. sie { 
B55 3k = P| (Xm Xm Dix — %j XE) 


Where B is the mass of the displaced fluid, 5,, is the Kronecker delta, and 
x, is the k-component of the centroid of volume in the body. 


3. MATHEMATICAL PRELIMINARIES 


In this section will be collected various mathematical theorems and 
results which will be required in the subsequent sections. 
3.1. Potential functions 

Let ¢ be a potential function which satisfies Laplace’s equation at points 
where no singularities are present, and which, in regions where there is a 
source distribution of strength a, satisfies Poisson’s equation 


= —4n0, (8) 


In the neighbourhood of a point source of strength m at the point with 
coordinates x,,, the potential may be expressed in the form 


¢ = ¢'+ “ ’ 7 ae (X;,— Xis)(%j— Xis)> (9) 


or, since ¢’ is analytic in the neighbourhood of x,,, 


m : od’ 
d= — +(¢’), +(4,-%,)( =) 2 
r. Ox;]; 


—“s 
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and similarly the velocity field is expressible in the form 
mn ; 0v; 
v= — —+(v’),+(x;-% AS) aay (11) 
rs s 
Qi v= -—-—. 12 
— Ss (12) 
In the neighbourhood of a point doublet of vector strength y at the point 
x,,, the potential may be expressed in the form 








. a Xiqg—% ° 
¢=¢'- : = ’ a = y mq = (x; ai Xia l(% — Xia)s (13) 
d d 
or, since ¢’ is analytic in the neighbourhood of x,,, 
bn : Cd’ 

d = — “4 ‘o (d da t (x; = wd(s-) Sore (14) 

‘7 OX; Ja 
_ ov, - 
v= —(3u,n,n,—pb,)+(%)at(%-—Xa(— }] + (15) 

ry 7 f : OX; 7 


3.2. Gauss’s and Green’s theorems 

Let 4(x,, Xs, v3) be a function analytic in a region R and on its boundaries, 
where the region RX is bounded externally by a closed surface S and internally 
by a set of spheres whose surfaces are collectively designated by S’. ‘The 
sense of the normals to the boundaries outward from the region R is 
taken as positive. 

We can now state Gauss’s theorem in the form 

| dn, dS = | - dr —| on, dS’, (16) 


where the integral on the left extends over the surface of the outer boundary, 
the last integral over the surface of the spheres, and the first integral on the 
right is a volume integral over the region R. Also, if (x1, %2, x3) is another 
function analytic in R and on its boundaries, we can state the second of 
Green’s theorems in the form 


i ae) A a7 a9 \ Z A ~ 
CW ( d . onMw o*dh ( us od ai 
(4 = —uJ =) dS = (eo: = — os = = dr— | o— — Y— dS ‘ 
\" On Cn, J\ 0X; 0X, ON; ON; . on on 


i 








or, if d satisfies Poisson’s equation and & Lapiace’s equation, 


Cus Ch : ; Cus od oe af 
(6 ah < | dS =4r | do dr—- (os —p=—]} dS’. (17) 


/ 
3.3 Boundary conditions and the Bernoulli equation 
Various partial derivatives of the quantities V’, and W defined in (4) and 
(5) will be required. These are given in this section. 
We readily obtain 
—— = — — +6€,0;, (18 
COX: CX, an ) 


If d satisfies Poisson’s equation, we have, when k = 1, 


= = —4no. (19) 
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Next, if ¢ satisfies Poisson’s equation, noting that dv,/éx, = dv,/cx,, and 
applying (19), we obtain 
ow C : 
s— = a—(-07,V + €,;;, dw) — 4107, . (20) 
CX, OX, 
Then, from (27) we obtain 


. , 
Eijk x; 7 ie € al = ae (x; Up y i+ J j U, —470x; On — Ekin x; Uv; oy | 
( . 


l 
But, substituting for V,, and simplifying, we have 


A 
ra Cc 
A 


Ei J j VU; = Ekim x; v, w») = €,;;(U, U;, + Exim W; ax, 


4. ADDED MASSES 
If we consider the kinetic energy T of the fluid due to the motion of the 
body when all the other boundaries and flow-producing mechanisms are at 
rest, we have the well-known formula 


2T = —-p | i dS. (22) 
From (1) and (2) it follows that 
2T = Aug Uy Ups (23) 


in which the added masses 4,9 are given by 


A 2 p ( d, Ng dS, A xB = A by. (24) 


ee) 


Can the added masses be expressed simply in terms of the singularities 
inside the body under consideration as some of them (A,,) were by ‘Taylor, 
Birkhoff, and Landweber for a body moving in unbounded fluid? ‘This is 
the chief concern of this section. In the following, two-dimensional and 
three-dimensional flows will be discussed separately. ‘l'wo dimensional 
flows are discussed not only because the use of function-theoretic methods 
enables one to solve the problem for two dimensions from an entirely new 
approach and in a singularly simple manner, but also because simple relation- 
ships between 4, and the singularities can be shown to be non-existent for 
pure rotations. 


4.1. Two-dimensional flows 

The complex variable z is defined, as usual, to be x,+7x,. With every 
velocity potential ¢, (x= 1, 2, 6 throughout the discussion for two- 
dimensional flows) one may associate a stream function %, such that the 
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complex potential is w,=¢,+i%,. Now, if s denotes the curvilinear 
distance along the body, we have n, ds = dx, , n. ds = — dx, , and since 
4) 14, 
cw ¢ om 
oe — = —n,, (25) 
cs on 


A. +tAys =p b fh, n, ds+ip fh, Nz ds = Ip d, dz 


. 


= -—Ip | ¢ w,dz—1t p wb, ds | = —Ip p w, dz—p ) sn, ds. 


. . 


Thus, 


Ay t+ By tt(Ayet+ By) = —ip p w, dz, (26) 


. 


in which 


B.. =p > x, n, ds, B=) ) von, ds. (27) 


By Gauss’s theorem, we have 


Bj; = Bs ;;, Be = — Bx, Beg= BX, (28) 


where B is the mass of the displaced fluid per unit length of the body, and 
v;, ¥, are the coordinates of the centroid of its area of section. 

Equation (26) readily yields the generalized ‘Taylor theorem. If there 
is a doublet of strength x, inside S, where in general j, is a complex number, 
there is a term y,/(s—%,) in #, and its contribution to the integral in (26) is 
27,1, by Cauchy’s theorem; and if there is a source of strength m, at 2,, 
inside S, there is a term —m,log(s—-,) in w,, and its contribution to the 
same integral is 27m, 2,7. ‘Thus, extending the result to distributed sources 
and sinks of strength o, per unit area of section, we have 


At By +t(Ayet+ By) = 2p] | 0, zs dA+X (m, 3, +1) (29) 


in which dA is an element of the area of section over which the integral 
extends. Equation (28) gives precisely the generalized Taylor theorem for 
two-dimensional flows. 

What, then, about the added mass for pure rotation? ‘The added mass 


in question is 


. 


fag =p ‘) d,n, ds =p .) ,(X| Mg—XyN,) ds = —p hy(x, dx, +X, dx5) 


= — p4#: p Ww, s*dz>+p ‘) W(X. dx,— x, dx3), 
where 3* is the complex conjugate of z. But, putting 


F(s) = | 2%”, ds, Fs) = | Xomeds, 
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integrating by parts, and applying (25), we have 


p b6(X_ dx,— x, dx) = - hg(X Ny + Xy Ny) ds 


- -> (1 Mg — Xqmy)(F + Fy) ds—f6(0) > (xy My + Xy My) ds 


= Af(xi + x3), — 24,(0)] p (x, Ny +X Nz) ds 


iy 
e 


—4 p (x7 +.x5)(x1 M+ Xy Ny) ds, 


where the index 0 denotes initial values along the path of integration. 
Applying Gauss’s theorem, we obtain 

i) pg(X_ dx,—x, dxy) = A[(x}+ x3)y—2,(0)] -—2 | (x7 +45) dA. 
We may take (0) = 0. Also, let r, be the radius of gyration of the area, 
and let 75 = (x7 +.%3)). Then we have, finally, 


A gg + Bog = pA(r,— r=) — pF: > We z%dz >. (30) 


It is now seen, by expressing w, in terms of its singularities, that (30) 
gives a linear relation between Age and their strength. If, for instance, 
among the singularities of w,, there are doublets of strength yu, at s = s,, 
each would give rise to a term 


. ~* 
— pp: > a dz P 


» 
~~ <; 





which, however, cannot be expressed in the form p,f(2,, 22, ...) independent 
of the shape of the profile; for otherwise, as a consequence of Morera’s 
theorem of complex variables, the integrand would be an analytic function, 
which it clearly is not. In this sense a simple relationship like those of (29) 
does not exist for A,,. The writers have laboured much and in vain in 
their endeavour to establish the missing simple relationships for pure 
rotations, and it was not until equation (30) was reached that they recognized 
the futility of their efforts. 


4.2. Three-dimensional flows 
For the fluid inside S and outside the small surface S’ around the interior 
singularities, we can apply (17) to obtain 


C 


Y i 
— d§ 
on 





A, = A XJ 


=p | dy n; dS = p | Dy 


_ dS +47p | x, 0, dr+p lore —¢o,n;) dS’, (31) 


Oo 


A 


« 








where o, is the distributed source strength corresponding to ¢,._ Here, as 
shown by Landweber (1956), 


<2" ; , 
p|x;=- dS = —p | x,n, dS = —B,,;. (32) 
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Moreover, for the source, 
| 2 — dS! = 4nEim, x, (33) 


For the doublet 1, , with components 2, , the potential is given by (13), from 
which, writing 
" Ody Og, 2x54 Mai Nj a Hei n; n; 


eo 3 2 
on ra T7G 








on : 


and noting that { x,(¢¢’/dn) dS’ vanishes in the limit as r, approaches zero, 
and fn,rj7°dS’ and fn,n,;r,? dS’ vanish because of anti-symmetry except 
when 7 = 7 in the last integral, we have 





- ad . (iM: 
| as ds = 2a , dS’ = - Hai» (34) 
J Cc : ~ 


where j is not summed in the second omiba of (34). Next, considering 
—f ¢,n, dS’, the contribution to it by a source is zero, since its potential 
varies inversely as the radial distance r, from it, whereas dS varies as r°. 
The contribution of a doublet is 


$’ “ o nN; , An a 

= \(¥- an n, aS = Uy; | a dS = 3% Pas: (35) 

Substituting (32) to 35); into (31), we obtain the generalized Taylor theorem 
A, + By; = 47p[fo. x; dr + (m,, x;,+ y;)]- (36) 


4.3. Evaluation of ‘ missing’ added-mass relations for an elongated body 
In this section simple expression for the added-mass coefficients for 
pure rotation in terms of the singularities within the body, will be derived 
for ellipsoids, and approximate ones for elongated bodies. 
From (17) we have 
i YY ( e ’ 
| d(x; ,.+x,n) dS = | $a (x, x,) dS 
1 0d 5, : ' ad 
= | 5 Xp dS +4r | ox; X;,, dT — | | 4, Ny +X}, Nj) —%; XK =| as" 


- Oo - ¥ a 
= | %; Hp ~ dS + ta | (ox, x, dr +X (mx; Xj, + psy Xp + by x) (37) 
; é : 
where the integrals over the spherical surfaces about the isolated singularities 
have been evaluated by a now familiar process and the subscripts s and d on 
the coordinates have been omitted in the last terms. Also for the first 
integral on the right in (37), we readily obtain from (3) and (16) the matrix 
of values 


pt ae es dS = B; 


3+t,0 


Bosiass ® ~Ple | (x°— x) dr, + (38) 
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where 2, j, k are ditferent and not summed. ‘Thus we may write 


| 


P| d(x a, +x%,2;) dS = By, . t4mpXx,, 


e (39) 
Dike | Oy %%, Ar+X (mM, X; Xp + fas XE + Pan Xj) 
\lso we have ; 
A, 3.i=? | py Nz..,;dS = pé; 55. | d, x, n; dS. (40) 


Consider the case of an ellipsoid rotating in an infinite fluid in which no 
other boundaries or singularities are present. ‘Then we have, from Lamb 
(1932, p. 154), when the coordinate axes are taken along the principal axes, 

Oy = CHa Ks; é,.- Cle e,, bg = Coa, Xs, (41) 
, C, are constants, and hence, from (39), (40), and (16), for the 
case x = 6, 1 = 3, we obtain 
Y a >/ PEE 
Cy Beg = Bag + 4 pXj06 


where C,, C; 


2 


6 
Ag ae Cs Bog ’ (42) 
or eliminating C,, 
Be ’ v 3 
“166° Rr r Boe a Oe t7pX196 . (+3) 
66 


Similar expressions may be written by symmetry for 4,, and A,,. The 
coefticients 4.,, 4,,, 4,,are zero for the chosen orientation of the coordinate 
axes. ‘hus we have found simple relations between the rotational added- 
mass coefficients and the singularities for ellipsoids, albeit not as simple as 
those in (36) 

When the x-axis of the ellipsoid is much greater than the others, and 
the x,- and x,-axes are nearly equal C, is very nearly 1, 44, <= —B,, from 
(42), and the left member of (43) may be written as 

> > 
Ass~ Aual “3 ! ) Big = Acct Biol a 1)- Big = Aca + Beg » 
66 66 
which is of the form that an extension of (36) would suggest. 

Next let us consider a body elongated in the direction of the x,-axis with 
nearly elliptical sections in the planes x, = constant. From (39) we obtain 
directly and exactly 


A B.. 4 2p | de Vo Ny dS = 42rpdX106 - (44) 


66 66 


For an elongated body the third term on the left is small compared with the 
other terms, so that only a small error would be introduced by assuming an 
approximate value for d, in that term. Since the section is nearly elliptical, 
let us assume ¢g== Cx, xy. Then 





a ee ee a ee — | xedr, 
- a Bas 
whence, substituting into (44) and simplifying, we obtain 
Bog 
i ets _— 2 
Aes Br + Bog = —4rpX jog. (45) 


66 
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The corresponding formulae for A,, and 4;; can be written down from 
symmetry. As in the case of the elongated ellipsoid, when the section is 
nearly circular we have also 

og + Bog = 4727p Xy26 (46) 
and similarly for 4;;. For A,, the term corresponding to the integral in 
(44) is of the same order of magnitude as the others, so that the errors in the 
foregoing approximations would be greater than for A;, or Ag, . 

Similarly we have, from (39), 


As,—Bigt 2p | bg X13 dS = 47pSaic , (47) 
in which the third term may be expressed approximately in the form 
2p | og x, Nz dS = 2pC | x? x. n3 dS = 0. 


Hence (47) becomes 
As, + Bsg = 47p Lsy¢- (48) 
An alternative relation for 4,, may be derived by interchanging the 
roles of the indices, and similar expressions may be obtained for 4,, and A,;. 
These as well as the relations for 4,,, 4;;, and 4,, may be summarized in 
the form 


B 


XX B’ 


daRr Ly LX -(A,,+B,,) = 7 7p sks (a = 3 a e- 4), 


+ (49) 
Ag+ Bg. = 47pe,;;, Ly, (RB =3+j, y=3+h), | 


in which 7, k are not summed and /, j, k are different. 


4.4 Linear and Angular Momentum 

It will be of interest to express the integrals p { du, dS and 
pe;; | ox;n,, dS in terms of the added-mass coefficients and the singularities 
within the body. By applying Gauss’s transformation to the region exterior 
to the body it is seen that the sums of such integrals over all the boundaries 
give the linear and angular momenta of the fluid. 

The velocity potential, in the form given by (4), may be further resolved 
by writing 

$o = b+ 9%; (50) 

where ¢/, is the potential of the part of ¢) due to external singularities, and 
¢, that due to internal singularities. We have then, from (24), 





p | dn, dS =u, Aytp | (45+$5)m, a5. (51) 

But, by Green’s reciprocal theorem and (17), 
| dn, dS = | x, 22 aS, (52) 
| ¢,, dS = | — dS +42 | o9x; dr+ i(, oi - om) dS’, (53) 


F.M. Y 
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where o, is the strength of the distributed singularities within the body 
corresponding to ¢y. ‘Thus, substituting (52) and (53) into (51), applying 
the boundary condition (3), and noting that the last integral over the spheres 
about the singularities is of the same form as that evaluated in (31), we obtain 


p | dn, dS =u, A.,t+47 | foo x, dr+X (my x; 4 10) (54) 


where m, and jo, are the strengths of the sources and doublet components 
within the body corresponding to ¢,. Hence, applying (36), we obtain 


p | on, dS+u, By, = tr) | ox, dr+X (mx “n)| (55) 


where o, m, and p, denote the totality of all the distributions and singularities 
within the body. ‘This shows that the value of the momentum integral 
depends simply upon the internal singularities. 

Next let us consider the angular momentum integral. Recalling the 
definition n,.; = €;,,%,n,, we have, putting 8 = 3+7, and applying (24) 
and (3), 


pe. | Ox,n, dS =p| (u,b,+b9)ng dS 


i , ., oh : - 
= U A.g—p | (J) +45) <5 dS. (56) 


X on 


But, from (17), we have 


: od ] 3 od), . ’ : od, , od vy 
| dg, "dS = | dg—~dS—4z | 034, dr+ (43 — — $5 $2) as’. 
On 2 , i on on 
Also, putting 3 = ¢,+¢3, where 4; is the part of dg due to external 
singularities, and 43 that due to the internal ones, we have from (17) and 
Green’s reciprocal theorem 


i 2 ~ to - i Vd ” Ad \ 
n OP 2 . CP , , OP n SPB ” 
| db, — dS | db > —— as + 4a | Oy db dr _ | db, — =n dS Fe 
on cn : : " On cn 
Ch " se ch ) Z 
bd,~— dS | d, — aS, 
on F on 


whence, applying the boundary condition (3), evaluating the integrals over 
S’ by the usual procedure, and substituting into (56), we obtain 


pe jl 


| px; ni; aS = u, As + tpl | (oa 8% dh.) dr + 
+X (mz db, — mo $3 — 63; Vo) + Ho vi) | (57) 


or, substituting for 4. from (36), and applying (7), 


pez | Oxyn, dS = Be, uj; X,+47p« | [o3(6) +u; x;) — 9 63] d+ 


+X [mg(d, + u; x;)— Mo >; — 143,(%o, — 4) +o; @3;] > +0; Agaigsj- (57a) 
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It will usually be convenient to choose the origin of coordinates so that 
the term #e;,u; x, vanishes, as is done in the equations of rigid body 
dynamics, where this term also occurs. 

For the elongated bodies considered in the previous section we can 
substitute the approximate values given in (49)for 4;_;,.;. Thus, applying 
the simpler of the approximate relations for 4g, we obtain 

pes, | bx, my, dS = Bes. u;X;,—w; Bs; 6+ 47p[| (a5 26 — 9% « ,) dr+ 
‘i eQ, ; 

+X (mg Q— my 3+ Me, a, + Boi Us )], (57b) 
= by tj Xj— Wy Ng Xz + Wy Xz Xy— Wz Ny Xo, 
and a similar expression for 7 = 2, with 

QO, = bo tu; Xj + wy Xe Xz — Wy Xz Xy— Wy Xy Xo. 

For the case 7 = 1, the simpler approximation for 4,, cannot be used so that 
the resulting expression for the moment of momentum integral would be 
more complex in form. 


"THE GENERALIZED LAGALLY THEOREM 


5.1. Force on the body 
From (5), the force on the body is given by 


F.=- | pn,dS = —p | ($ _ Wn dS 


1; : 
p33 dn, dS +p | Wn, dS, (58) 


where the prime in ¢’ denotes that the variation with time is relative to a 
moving coordinate system. But, from (16) and (20), 


OW 





| Wn, dS = | <— dr— | Wn, ds’ 


al 3 ig V, +4, beng) de~4er | om, dr—| Wn, ds’. 
Also, applying (4) and (16), we have 
| = Vit ein box) dt = €ig, OK | gon; dS + |e V+ <€ jn poox)n, dS’. 
iaais: substituting into (58), we obtain 


F, = -»= | on, dS + pe, 7. w;. | gn; dS —4zp | ov; dr+ 
+p|[(—2; V,+€;;,¢0,)n;— Wn] dS’. (59) 


The sum of the first two terms of (59) is seen to be the absolute time 
derivative of the momentum vector integral — p J én,dS for which values were 


Y2 
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given in (54) and (55). Thus, evaluating the integrals over S’ by the usual 
, 


procedure, we obtain the expression for the force 


‘ d d . e | 
if dt \s Iai) ‘ap Al Oy X; dr +X (my x; + Hw) | Z 


tr] | Ov ; dz + > (» rom TOLL, t pj =) | (00) 


or, if (55) is used in (59), and (7) ‘s applied to evaluate 





‘ 1 W;X,)= Bu, , 


9 » » » 
u,B 7 u; Bj, +; B: - Blu ,;+ 


where u, is the -component of the velocity of the centroid, then 


, du d a ) 
. dt ae a | ox, dr + X(mx, -u) | 


4p] | ov. dr- ( v. —snou,+4 Ms! ‘ 61 
| l > (™ TOL, wi) | (61) 


If it is desired to use relative rather than absolute time derivatives, we may 
employ alternative expressions of which a typical one ts 
d 


dt om) 


dt (ox,) + O€ ; 57, W Vp. 
Furthermore, it should be noted that the strengths of the singularities 
occurring in (60) and (61) are those due to the superimposed effects of all 
the velocity components of both the body and external boundaries and flow- 
producing mechanisms, in contrast with the singularity strengths corre- 
sponding to unit magnitude of a single velocity component of the body 
which occur in the generalized Taylor formulas for the added masses. 

It isimportant to observe that, in computing 7; and cv; Cx, in (60) and (61), 
the contributions from all the internal singularities may be omitted. The 
reason for this is that the mutual contributions of the velocity fields of a 
pair of internal singularities to the expression for the force on the body are 
equal and opposite and hence annul each other. ‘This introduces a signi- 
ficant simplification in the calculation of the force from these equations. 

Equation (61) is essentially equivalent to the corresponding result by 


Cummins (1953). It ditfers from it in the following respects. 


(1) A new term, 1l67%pcu, 3, has appeared. This did not occur in the 
treatments of Lagally or Cummins because they did not consider the 
case of simultaneous occurrence of distributed and isolated singu- 
larities. 

(2) The inertia term in Cummins’ equation (8) and the terms of his 
equation (48) for the vector F; have been replaced by the first term 
of (61), which is seen to represent the inertia of the displaced fluid. 


(3) The expression for the force has been expressed in a much more 
compact form. 














1 
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5.2. Moment on the Body 
The moment on the body is given, from (5), by 


M, = — x | px; n, dS = —pe;;,. It; - W)x n, dS 


( 


1 ° 
— pe | 7 bx,n, dS — |Wx,n, as |. (62) 
dt’. ; 
From (16), (21), and (4), we obtain 
. oe 
€, | Wx,n, dS = «, | RY = dr—|Wx,n, as’ | 
= —€ fu | dn, dS+w, | dny.;, dS+4r | ox, 0, dr |= 
B 
en | (xj vy Vin, +u,; dn, +e, dns. ,+ Wx,n,) dS’. 


Hence, evaluating the integrals over S’ by the already frequently applied 
procedure and putting 


i; e* 
Ennl sp | ox, dS+w, | dns., dS | = Poe | dxjn, dS, 
‘Lat ; “— 


the expression for the moment becomes 
dt . 
where 1/,,(7) denotes the Lagally moment for steady flow, 


; Ov, 4 
M,,(v) = —4pe, | ox; t, at +> (ms Up t py UgtX; Mpa —ATON ux) | 
CX, 


This ditfers from the expressions derived by Lagally and Cummins for 


M, = —pe iE ! bx, ny, dAS+u, | dn, as | M,,(v), (63) 


steady How in the appearance of the last term, because they did not consider 
the case in which distributed sources and isolated doublets are simul- 
taneously present. Except for this term, (63) may be shown to be equi- 
valent to the result derived by Cummins for unsteady flow. Cummins, 
however, did not succeed in expressing the moment of momentum integral 
in (63) in terms of singularities, so that its time derivative occurs explicitly 


in his final result. As in the case of the Lagally force, the contributions to 
from internal singularities need not be considered in computing the 
Lagally moment since they annul each other in summation. 


Now, substituting for the momentum and moment of momentum 
integrals in (63) from (55) and (57a), and noting from (7) that 
d _ _ at, 
Be... — (4; %,) +€,; 4,4, BL, = Be. x : 
at dt 
ve obtait 
u a d 
VW Be... x - 1 ) [oa(d, +u,x;)—o, 63] dr+ 


Y [m3(d,+u x;)— my b3 — 43,(%;— Uj) + Mo; aj] + Miz(v—u), (64) 


a 
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where 8 = 3+172, and M,,(v—u) denotes the value of the Lagally moment 

when the velocity at the singularity relative to that of the body is used. 
Sor elongated bodies with nearly ellipsoidal sections, substituting for 
‘rom (49), we can obtain an approximate but complete expression 


for )\ 1ent in terms of singularities. Thus, applying (57b), we have 
du, a d{; 
M, = Bea, X,—« + —(w,; Bs; e)— 4p —% | (og Qe— 99 4,) dr 4 
3 JERip dt ap 3+i, ;) Pat ik 6 76 0 Pe) 
. 30, | 
+i | m, OQ, — my $y + Mea +o vi | -+M,(v—u). (64a) 
Cx ia: . 


A similar expression for ..W, can be written down by symmetry. In the 
application of (64) and (64a) it will usually be convenient to make the first 
term vanish either by choosing the origin at a point of zero acceleration, or 
at the centroid, or at a point whose acceleration vector passes through the 


centroid. 


6. COMPARISON OF THE RESULTS OF KIRCHHOFF, TAYLOR AND LAGALLY 

Kirchhoft’s equations of motion of a body through a fluid express the 
force and moment acting on a body in terms of its added masses. Since 
the generalized Taylor theorem expresses the added masses in terms of the 
singularities, it appears that, by substituting these expressions for the added 
masses into Kirchhoff’s equations, it might be possible to derive identical 
formulas to those derived above for the forces and moments in terms of the 
singularities. 

First let us consider one of Kirchhoff’s equations (Lamb 1932, p. 168), 
for the case when the fluid is disturbed only by the motion of the body, 

; d oT oT oT , 
= ge ee ee Bs (09) 
where 7 is the kinetic energy of the fluid and the prime denotes that the 
time derivative is taken relative to a moving coordinate system. ‘Then 
F, ee u,(u; A,3—Ug Ay) = — © (ta Aaa) (66) 
On examining the expression for the corresponding term in (60), it is seen 
that the terms in the first bracket vanish since, in the present case, there are 
no external singularities and consequently no internal images of them. ‘The 
terms in the second bracket vanish because the mutual contributions of the 
internal sources and doublets to the force on the body annul each other. 
Consequently, (60) and (66) are seen to be in agreement. 

More generally, if it is supposed that other boundaries are present, but 
that only the given body is in motion, the force on the body due to the fluid 
may be obtained by applying Lagrange’s equations (Lamb 1932, p. 188) in 
the form 


ae 2T = A, gu, Ug, (67) 
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where €, (¢=1,2,3) are the coordinates at the origin of the coordinate system 
attached to the body, relative to an inertial system. ‘Then 


d 
dt 


0A 


. l 
F, = (u,, Agi)+ 5 i? Uy Ug. (63) 





Comparison of this result with (60), in which the first bracket vanishes in 
the present case since the external boundaries are at rest, shows that the 
two expressions for the force would be identical in form only if 


oA ; ov, 
ws = fee =" i 
=r Uy Ug = — Spf | ov, dr+ > (om i 37 OH, + Mya ) , 


CC 





i 


or, putting o = u, o,, U; = Ug V,, etc., only if 


cA ; . ( Ue\ 
xp ) = = = hes. g 
“3 = —$7 | a, Ug; dt + > (m, Upi — 37M yi 9B + My oe )] (69) 








€ 


But substitution of the generalized ‘Taylor formula for the added masses in 
terms of singularities into (68) does not seem to yield identically the genera- 
lized Lagally formula for the force; in fact, such a complete substitution 
could not be made since not all the added masses can be so expressed. ‘Thus 
the method of Kirchhoff-Lagrange appears in general to furnish an alter- 
native (and more limited) method of computing the forces on a body. 
Nevertheless, (69) must be valid, since the force may be obtained by either 
method, so that we have expressions for the gradients of the added masses 
in terms of the singularities. 

It will be instructive to illustrate both methods by computing the force 
on a sphere A of radius a moving with velocity u, along the line of centres 
away from a fixed sphere B of radius 6. Let c be the distance between 
centres at a given instant. Using the method of successive images, we 
begin with a doublet of strength py, = $u, a* at the centre of A. Its first 
image in B is a doublet of strength 

b3 u, a®b3 
P= F 73 Ho a sis 33 


at a distance = c—h?/c from the centre of 4. ‘This gives a second 


Hotbe _ ; 3.a*b* 
A,, = —B,,+4p a = frail 1 + (2-8 ’ 
and (68) gives for the force 
4A, d 67puz a®b3c 


d 
F,=- Frau Ay)+ 2%; =e 7 (ei An)— “(2—hy * 
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In order to apply (60) we also need the velocity v’, due to the first doublet 
image p1,, along the line of centres, 
v = 2u,(x — b/c), 
where x is distance measured from the centre of B. Then 
ev’ — 3a, abc 
ox (xc—B®)i° 





In order to obtain the same order of approximation as above, we need com- 
pute only the force on the original doublet. We obtain from (60) 


; d O7 d 6zpuiz a®bec 
F,=- ay Ay) —4npta ) aie ets Fre Aj,)- (ey 


“ , 
Jey 
— 
Ox 


which agrees with the result above. 


‘The present work was conducted at the lowa Institute of Hydraulic 
Research, State University of Iowa, under contract Nonr 1611 with the 
Office of Naval Research. ‘The authors wish to thank Milton Martin and 
Manfred Bottacini of the Institute for their careful reviews. We are also 
pleased, here, to mention George Weinblum of the University of Hamburg, 
that most inspiring teacher, who pointed out the power of the Lagally 
theorem and new fields of research for many of us. 
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SUMMARY 
Ward’s slender-body theory of supersonic flow is applied to 
bodies terminating in either (i) a single trailing edge at right 
angles to the oncoming supersonic stream, or (ii) two trailing edges 
at right angles to one another as well as to the oncoming stream, 
or (iii) a cylindrical section with two or four identical fins equally 
spaced round it. ‘The wave drag at zero lift, D, is given by the 





expression 
D 1 Lee he Vit vi ae? 
UE * i | ; | : log 2 S"(s)S"(2) dsdz — 
SL) (t l S’(1) l 





(M2 — 11m + */> 

where / is the length of the body, 6.the semi-span of the trailing 
edge (or length of trailing edge of a single fin), and S(z) is the 
cross-sectional area of the body at a distance z behind the apex. 
The constant k depends on the distribution of trailing-edge angle 
along the span for each trailing-edge configuration. In case (i) it 
is 1-5 for a uniform distribution of trailing-edge angle and 1-64 for 
an elliptic distribution. In case (ii) it is 1-28 for a uniform distri- 
bution and 1-44 for an elliptic distribution. Study of case (iii) 
indicates that interference effects due to the presence of the body 
reduce the drag of the fins. For example, with a uniform distri- 
bution of trailing-edge angle, k for two fins falls from 1-5 in the 
absence of a body to 1:06 when the body radius equals the trailing- 
edge semi-span, while k for four fins falls from 1-28 to 0-45 under 


| | log | Sana S"(z) dz T 


= < log 


the same conditions. 

Where ordinary finite-wing theecry is applicable, the present 
method must agree with it for small (M?—1)"%b//, and this is 
confirmed by two examples (§ 3), but within the limit imposed by 
slenderness the present method is of course more widely applicable, 
as well as simpler, than finite-wing theory. 

It is not known experimentally whether slender-body theory 
gives accurate predictions of drag at zero lift, for the shapes here 
discussed, under the conditions for which on theoretical grounds 
it might be expected to do so. It should be noted that, although 
tests have not vet been made on ideally suitable bodies, no clear 


w 
wn 
“J 
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indication of the variation with log(M?-—1) indicated by the 
formula has ever been found. 


1. INTRODUCTION 
Ward’s ‘slender-body’ theory (1949) applies to supersonic flow about 
‘bodies’ (here understood to include wings and wing-body combinations) 
whose surface is headed by a pointed apex and makes everywhere an angle 
with the undisturbed stream small compared with the Mach angle. Then 
the right-hand side in the linearized equation of motion for ¢, the disturbance 
potential*, 





= (M2-1)— (1) 


RoR ae 
(in cylindrical coordinates with the z-axis parallel to the undisturbed 
stream), may be expected to be small near the body compared with the 
left-hand side, provided that the shape of the body varies smoothly with z. 

Accordingly, ¢ may be determined in this region as a solution of 
Laplace’s equation, in a plane perpendicular to the undisturbed stream, 
which satisfies the boundary condition at the body surface. ‘The solution 
that satisfies 


; “Sez) l M?—1)2 _ i 
D ae logr + aan log at! a8 S (s) ds + Of ) (2) 
Lol 271 0 Az Ss) Y 
is r- >» © must be selected if it is to join smoothly, far from the body, on to 


a solution of (1) which is zero upstream of the Mach cone from the apex. 

This paper is concerned with bodies which terminate in a straight 
trailing edge at right angles to the stream, with a small but finite trailing- 
edge angle (which may vary along it). ‘This class of bodies includes delta 
wings. Bodies terminating in two such trailing edges at right angles to 
one another, or in a cylindrical section with two or four fins attached, are also 
treated. 

Ward and others have used his theory to obtain the lift, moment and 
drag due to lift for such bodies. In this paper their drag at zero lift is 
investigated. ‘lhe Ward theory is applied without change, so it may be 
asked why a separate paper is necessary. This is partly because Ward has 
suggested that the theory is inapplicable to. cross-sections with large 
curvature, and the trailing-edge cross-sections described above have infinite 
curvature at the tips One might hope, however, that this suggestion of 
Ward’s would be correct only if fluid flows round the region of high curvature 
(leading to specially low pressures there), and that in cases like those con- 
sidered here, where this is not so, the theory can still be alppied. 

Another reason is that many people misinterpret Ward’s theory and 
suggest that it equates drag at zero lift to that of a body of revolution with 
the same distribution S(z) of cross-sectional area with distance from the 
apex. For the bodies here considered, this ‘equivalent body of revolution’ 

* Such that U(2+4) is the velocity potential, where U is the velocity of the 


undisturbed stream. 
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would have a rounded rear end, which might lead one to suppose the theory 
inapplicable, especially as the formal expression for drag calculated for this 
body shape would be infinite. 

However, Ward’s actual drag formula for a body of length / is 

D 


ipl? 


1 rl 
=a-}| | log- x S"(s)S"(z) dsdz — 


0 ie —~ 6] 


Ps ‘l eG 
ii si) | log i S"(2)ds-($¢ >, ") (3) 
i Ov 


» = ~ 
27 l z=l 








~ 0 - 
(Ward 1949, equation (37); Ward 1955, equation (9.8.7); here, dv and dr 
are line elements in the plane of a cross-section, normal and tangential to 
it respectively; the ‘base pressure’ term has been omitted). Now, the 
first two terms in (3) are independent of the shapes of cross-sections, but the 
last term, an integral round the rear cross-section, is not. It is because this 
term vanishes in so many cases at zero lift that the theory of the equivalent 
body of revolution has been assumed to be true more generally than it really is. 
This third term is calculated below for bodies terminating in a straight 
trailing edge, leading to the drag formula* 





D Lf; 
a | | log = S”(s)\S"(s) dsdz — 
»pU~ <7 0. 0 S—8 
S(D | Dee ee l : 
on = R log * aap S (x) dz ss —" log (M2— 1)" +k >, (4) 


where 26 is the length of the trailing edge and the constant k is given by the 


| y)dy) ( 


} 


equation 


- } 


~o b 2h 
k= | | loo ———— ef vye(Y)dyvdY ( | 
vy a : 
where e() is the trailing-edge angle at a distance vy from the middle. For a 
uniform distribution of trailing-edge angle k = 1-5, while for an elliptic 
distribution k = 1-64. 

Equation (+) is directly applicable to delta wings when the parameter 
(.\/?—1)'*b7 is small. Experimental work by Love (1949) and others 
indicates that this is perhaps the only region in which any linear theory is 


wi 
— 


adequate for predicting the wave drag at zero liftt, although some results 
f Herbert (1951) discussed in § 3, where a particular wing shape introduced 
by Squire (1951) is considered as an example, are a reminder that it should 
not be applied with V/ too near unity. 

In $4 a similar analysis is applied to other trailing-edge configurations, 
leading always to the same form (4) for the drag with different expressions 
for k. The extent to which a ‘dart’ formed of two slender delta wings 
it right angles (with a common line of symmetry) has a drag approximately 
equal to four times that of each separately, as suggested by the theory of the 


* The replacement of the 27 in the second term in (3) by z in (4) ts not a misprint. 
+ Perhaps because outside this region either doubtful assumptions about ‘ leading- 


edge forces’ have to be made, or else the flow norma! to some edge is nearly sonic. 
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equivalent body of revolution, is discussed; the true slender-body theory 
drag is somewhat less than this. 

An analysis (§4) of bodies terminating in cylindrical sections with 
slender fins having trailing edges at right angles to the stream shows that 
interference effects due to the presence of the body tend to reduce the drag 
of the fins. It had been widely supposed that this would not be true in 
the case of two fins, but the present result, which is subject to the usual 
limitations on accuracy of slender-body theory but can hardly be without 
qualitative validity, indicates that caution is necessary in inferring the drag 
of delta wings from flight tests on cylindrical projectiles which carry them 
in the form of fins. 


2. GENERAL THEORY FOR A STRAIGHT TRAILING EDGE PERPENDICULAR TO THE 
UNDISTURBED STREAM 
In the plane, at right angles to the stream, which includes the trailing 
edge, we choose Cartesian coordinates (x,y) such that the trailing edge is 
x= 0, —b <—y <b. The boundary condition on ¢ in this plane is 
ad _{— dey) (x= + 0) 
Ox +4e(y) (x = — 0) 
where ¢(y) is the trailing edge angle at a distance y from the centre. 
A solution of Laplace’s equation in the (x,y) plane satisfying this bound- 
ary condition is obtained by a sink distribution, the sink strength in the 


(-b<y <b), (6) 


interval dy being e(v) dy; this gives 


“bh & Y ; ' = 
des -) Ble t ety Peery ic. (7) 
h “77 
where C is an arbitrary constant. ‘Then as r = (x? + ?)!?_- ow, 
e(Y Sd 
d+ (| a logr b OO ogr oP (S) 
since evidently S’(/) | e())d¥; and hence, by (2), 
l (Mz—1)12 
( ot ge XI 2) Ps (x) dz. (9) 


Now, Ward's integral, which has to be taken around the boundary of 
the rear cross-section, is 


{ Co bd b ch (* C db \ 
~ d— dar | db = dy+ | o— d\ 
: ol } OX o h ox ) 


| [° e(y)e(Y)logly— Y|dydY, (10) 


at Be } 


CS'(1) 4 


by (7). A combination of equations (3), (9) and (10) gives at once the drag 
formula (4), with the expression (5) for the constant k. This expression 
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for k has the merit that it is independent of b or of the scale of e(y), depending 
only on the distribution of trailing-edge angle along the span. 
Thus, for a uniform distribution (¢(7) = constant), 
‘i 


1/1 2 
sme i | | es ¥ =¥) dydY = 2 = 155, (11) 
See ee ee 


while, for an elliptic distribution, 
‘1 Z 
| (-yfed-Y lee yo 


= log ++} = 1-636. (12) 


EXAMPLES : APPLICATION TO DELTA WINGS DESCRIBED BY SQUIRE AND 
BY PUCKETT 
Squire (1951) considers a delta wing all of whose cross-sections by 
planes normal to the stream are ellipses. ‘he major semi-axis is bz// and 
the minor semi-axis is 272(1—://), where 7 is the thickness-chord ratio of 
the central section of the wing*. Slender-body theory is applicable to 
such wings if (M?—1)!*6// is small. 
The cross-sectional area distribution for this wing is 
S(z) = 2(bs/l) (2r2(1—2/l)}. (13) 
The trailing-edge angle is equal to 47 at the centre and varies elliptically 
over the length 24 of the trailing edge. Note also that 


b 32\ 
S'(l) = -—2n7b, S"(s) = 4u7-(1- — }. 
(!) $a) = dor (1- F) (14) 
Substitution of these values in (5) gives for the three terms (respectively) 
D ‘ i l 
——. = 2 72h? — 10 72h? + 20725? ———_—____. + 1-636 } 5 
pO? 3 10277 Tb | 18 Th sand (15) 


the value of k appropriate to an elliptic distribution of trailing-edge angle 
being inserted. ‘The drag coefficient based on the area /b of the wing 
planform is therefore 

D 50 { 


l ) 
pee = Dr? 4 lo, 4 REED. 
Cp = co. *- mF log apy ian (16) 


This agrees with the asymptotic form of Squire’s expression for C,, for 
very small values of (M?—1)!%b//, but is somewhat smaller for moderate 
values, as figure 1 shows. (Note that the form of Squire’s expression 
depends on the precise hypothesis used to get over the notorious ‘leading- 
edge force’ difficulty, but that the asymptotic form for small (M?—1)!2b/1 
is independent of this, as no large pressure coefficient at the leading edge 
arises under this condition.) 


* This central section is biconvex, but there is a round-nosed aerofoil section at 
all other spanwise positions. 
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The free-flight measurements of Herbert (1951), in the range 
(M?—1)!2b/1 < 0-4 
where the slender-body theory can reasonably be applied, are well below 
the theoretical values, but the wing used by him was so far from slender 
(in fact b// = 1) that the condition mentioned holds only when V/ 1-08, 
so that the problem has become a transonic one. It is hoped that 
transonic version of slender-body theory may one day be developed*, but 


this has not yet been done. 





Figure 1. Drag results for Squire wing 


As another example, consider a delta wing with a double-wedge aerofoil 
section of uniform thickness-chord ratio + (Puckett 1946). Then the 
cross-sectional area distribution is 

; 2rbz?/1 z 1), “ 
S(s) = <5 a a (17) 
27b(4s —1—3s7/l) (zs > $1). 
The trailing-edge angle is 27 all along the span, and 
4751 (2 < Jl) 
S‘(/) 4752S (2) = 5 ar/+\ 18 
( 4) =) _49ebid (x > AD. aia 


Substitution in (5) gives for the three terms (respectively) 


) 872b2 /, 82h? 
= ; — (: t +log2) sills ane (2 + +log 2) + 

- 877b? | l P 19) 
+ ecu Tap 0) 


fact that equation (1), on to a solution of which the harmonic function ¢ must join 
f the ‘ transonic’ term 


This would require the modification of the second term in (2), to allow for the 


far from the body, is modified by addition o 
- 


(y-DIP = 


to its right-hand side. 
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the value of & for a uniform distribution of trailing-edge angle being 
inserted. ‘The drag coefficient is therefore 
877.5 l ) 

edad 8 (1) Tf? 
in agreement with the asymptotic form of Puckett’s expression for (p 
for very small values of (/?—1)'*b/l. The difference is a term 

O[(M?— 1)b?//?] 

inside the curly brackets in (20), and Ward’s arguments indicate that this 
will generally be the case. As with the Squire wing, Puckett’s expression 


(20) 














fe) f 02 0-3 04 0-5 0-6 07 08 0-9 0 
(m2-1)2 b/2 
Figure 2. Drag results for delta wing with double-wedge section. 


for the wave drag rises above expression (20) as the condition of slenderness 
is removed, as figure 2 shows. 


+. ‘THEORY FOR OTHER TRAILING-EDGE CONFIGURATIONS 

One interesting configuration is a body terminating in two trailing edges 
at right angles to one another as well as to the oncoming stream; for 
example, a projectile of revolution tapering to zero radius at the rear end 
and stabilized by cruciform fins. ‘This can be treated by methods similar 
to those of § 2. 

We suppose that the distribution of trailing-edge angle is the same on 
each trailing edge, and choose the x-axis and y-axis to lie along the two 
edges. Then there are two boundary conditions on ¢, one being specified 
by equation (6), and the other being 

ch _ {—$ex) (y= +0) 


ey - | Ss L(x) (y = —V) if > b<x- a b). (21) 
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A solution is obtained by a sink distribution as in § 2, giving 
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with C given by equation (9) as before. Note that now 
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where four double integrals have been combined into one by taking the vari- 
ables of integration as y and Y ineach. A combination of equations (3), (9) 
and (23) now gives the drag formula (4) with the new expression 
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fork. As before, k depends only on the distribution of trailing-edge angle, 
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not on its magnitude or on 4. It tends to be slightly less than its previous 
value (5) because on the average (vy? + }?)'? exceeds |y— Y 
Thus, for a uniform distribution of trailing-edge angle, 
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(compare 1-5 in the case of a singie trailing edge), while for an elliptic 
distribution 
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(compare 1-636 in the case of a single trailing edge). 

The fact that & is less than the corresponding value for a single trailing 
edge means that the drag of a pair of delta wings superimposed with their 
apexes and axes of symmetry coinciding, and their planes at right angles to 
one another, is not quite four times that of each separately. For example, 
a ‘dart’ of this kind made of two Squire wings has 
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where the drag coefficient is based on the sum of the wing planform areas, 
and the constant is different from that in (16) owing to the different value of 
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k which is appropriate. However, for very slender wings the logarithmic 
term is the most important in (16) and (27), and then the interference 
between the wings almost doubles the drag of each. 

We may consider also bodies which terminate in a cylindrical cross- 
section and a pair, or two pairs at right angles, of fin trailing edges. Let R 
be the radius of the cylindrical cross-section and 4 the length of each trailing 
edge, and suppose that the trailing-edge angle is «(qg) at a distance g =r—R 
from the cylinder. ‘Then a solution for ¢ in the plane, at right angles to 
the stream, which includes the trailing-edges, is obtainable by a sink 
distribution as before, but image sources and sinks in the cylinder have to 
be included as well, giving for the case of two fins placed opposite to one 
another 
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A combination of equations (3), (9) and (29) now gives the drag formula (4) 
with the new expression 
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fork. ‘This expression depends as before on the way in which the trailing- 
edge angle is distributed, as well as on the ratio R/b = a, say. 
For a uniform distribution of trailing-edge angle, 
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and this integral is straightforward to evaluate if the heute of the curly 
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bracket is expanded in series. ‘his gives values of k as a function of 

= R/b as in table 1 (see also figure 3). It is seen that the effect of 
interference between fins and body is to decrease the drag of the fins by 
decreasing k. In the limit of very large R/b (which is of more theoretical 
than practical interest) k becomes }—log2 = 0-807. ‘This may be thought 
surprising on the grounds that the drag of a wing composed of only the 
two fins (the case R/b = 0) ought to be reproduced for very large R/b, when 
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each fin, being accompanied by an image in the nearby plane boundary otf 
the cylinder, should have half the drag of such a complete wing. However, 
this argument neglects the effect on a fin of the image sources at the centre 
of the cylinder, and of the sinks (and image sinks) representing the other 
fin. It is the fact that the latter are twice as distant as the former, though of 
the same strength, which is responsible for the limiting result 
(R) ny — o@ = (R)pp m y— log 2. 

Comparable behaviour is to be expected for other distributions of 
trailing-edge angle, and the limiting result just quoted is easily seen to be 
quite general. 
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Table 1 


Note that the limiting process R/b—- « should be thought of as achieved 
by reducing 6 for given cylinder radius R and fin length. If instead R is 
increased indefinitely, keeping 4 and the fin length fixed, the two fins must 
cease to interfere with oné€ another after atime. At the same time, slender- 
body theory will clearly have ceased to be applicable. ‘Thus, there are 
certainly circumstances in which the interference here discussed does not 
take place. However, in the many practical cases when the two fins do 
interfere, the present theory, although only approximate, makes it appear 
likely that interference will at least lower the drag. 

Consider a numerical example. Herbert (1951), in the experiments 
quoted in figure 1, attached his delta wing in the form of two fins to a body 
with R/b = 0-16. ‘This, by table 1, should reduce k by 0-15, and hence 
reduce the drag coefficient C,, (for a Squire wing with 7 = 0-06 and 6// = 1) 
by 27(0-06)?(0-15) = 0-0033. This reduction is insufficient to remove the 
discrepancy noted in figure | (and there ascribed to transonic flow condi- 
tions), but such a relatively big change for such a moderate value of R/b 
indicates that the effect here discussed is by no means negligible. 

Passing to the case of a body with cruciform fins, we easily obtain by 
similar means the formula 
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For a uniform distribution of trailing-edge angle this is 
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where « = R/b as before. This expression gives values of k as function 
of R/b, as in table 2 (see also figure 3). 


\ _- Asymptotic value 


Rib ic 


Figure 3. Values of the constant k in the expression for wave drag at zero lift given 
»¢ slender-body theory, for a body terminating in a cylindrical portion of 
be slender-body theory, f body t t vlindrical portion of 
radius R, and two or four fins, the length of whose trailing edge is 6 in each 


case. The trailing-edge angle is taken as uniform 


Now, the behaviour of k exhibited in table 2 is crucially different from 
that in Table 1, for k decreases indefinitely as R/b is increased. ‘This was 
to be expected, on the grounds that with four fins conditions for very large 
R/b become fundamentally different from those for small R/b. Then, 
quite apart from the remoter interference effects, each fin by virtue of its 
image in the body would have half the drag of a complete wing made up of 
two fins. ‘The total drag would therefore be expected to be about twice 
the drag of such a wing, whereas for R/b = 0 slender-body theory predicts 
(as shown earlier) a drag only slightly less than four times the drag of a 


single wing. 
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Table 2 


Hence, the steady decrease of k as R/b increases, shown in table 2 and 
figure 3, represents as far as slender-body theory is capable the gradual 
tendency to a state in which the fins cease to interfere with one another and 
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the drag is therefore twice that of a wing made up of two of them. 


stages ol the process cannot be repre sented 


initial trend may well be indicated fairly ac 
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REVIEWS 


The Theory of Hydrodynamic Stability, by C. C. Lin. Cambridge 
University Press, 1955. 155 pp. 22s. 6d. or $4.25. 


Hydrodynamic stability is no newcomer to the literature of fluid 
mechanics. It was studied by many of the 19th century scientists, who 
recognized its important part in the neat structure we now call classical 
hydrodynamics. ‘The question of instability in vortex sheets was considered 
by Stokes and Helmholtz, and a number of other special problems were 
dealt with successfully. It is instructive to recall also the care and skill 
with which scientists of the last century examined the general question of 
the modes of disturbance possible in a continuous dynamical system. 
Intuitive concepts of stability are not easily extended to a continuous medium 
in motion, and a contemporary physicist regarding the present subject 
from a non-expert’s viewpoint might still attach foremost importance to 
this basic issue. Thus there may be some who would regret that the book 
under review says very little about it. "This is scarcely cause for complaint, 
however, for the book is admittedly preoccupied with the intricate mathe- 
matics of modern theories. 

In this as in so many other subjects, the first definite steps towards the 
modern science were made by Lord Rayleigh, who attacked the general 
problem of the stability of plane flows in which the velocity and vorticity 
have a continuous distribution. His work culminated in the well-known 
theorem which states that a flow is stable when its velocity profile possesses 
no point of inflexion. It is interesting that he commenced this work a 
few years before Reynolds made his important experiments first revealing 
transition to turbulence in a pipe, although this discovery naturally 
encouraged him in his later researches. (Strangely enough, a critical 
Reynolds number for Poiseuille flow has not been established theoretically 
even to this day.) ‘The problem of explaining the mechanism of transition 
is still the chief stimulus to research in hydrodynamic stability, and our 
answers to the crucial question how random turbulent fluctuations develop 
from oscillatory laminar flows remain no more than tentative. Neverthe- 
less, in those cases where the initial laminar flow is only slightly disturbed 
and the solid boundaries are smooth, an essential first stage in the process is 
without doubt an instability, which results in the amplification of oscillations 
whose wave numbers are contained within certain bands. 

The confliction of the inferences from Rayleigh’s theorem and from 
Reynold’s experiments led the earlier workers to the conclusion that 
viscosity plays an essential part in the development of instability; and, 
after several unsuccessful attempts to explain transition by reference to the 
energy of the disturbances, the basis of the modern theory was laid by 
Orr and Sommerfeld, who derived the equation governing the propagation 
of small wavy disturbances in otherwise steady laminar flow. Subsequent 
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work on the solution of this equation was principally directed, on the one 
hand, by G. |. Taylor to the flow between concentric rotating cylinders and, 
on the other, by Heisenberg to the much more difficult case of two-dimen- 
sional Poiseuille flow. ‘Taylor performed experiments which substantiated 
his theoretical results with extraordinary completeness; but Heisenberg’s 
analysis contained certain mathematical obscurities, and, besides the lack 
of any experimental support, his conclusions were contradicted by other 
investigators. ‘The extension, by Tollmien and Schlichting, of Heisenberg’s 
analytical methods to the case of boundary-layer stability carried with it the 
doubts attached to his earlier work, and was also subject to the criticism 
that the growth of the boundary layer had been neglected without wholly 
satisfactory justification. The controversy led to a decided schism 
between workers in Germany, where development of stability theory 
continued, and those in other countries who preferred an alternative theory 
proposed by ‘Taylor. The conflict was finally resolved in favour of 
stability theory—at least, in circumstances where the disturbance of the 
laminar regime is small—by the experiments of Schubauer and Skramstad. 
Their investigation revealed boundary-layer oscillations whose structure 
was in remarkably good agreement with the predictions of Schlichting, 
for which further confirmatory evidence was subsequently found by 
Liepmann. ‘The experimental vindication of stability theory eventually 
led Lin to reconsider the mathematical foundation of the subject, which he 
was able to clarify in large measure. Indeed, the present status of the 
theory owes a great deal to Lin’s efforts. 

There can be no reasonable doubt about the usefulness of Professor 
Lin’s monograph, for it is the first major survey of the theory of hydro- 


dynamic stability to appear in book form. For a complete account of the 
physical aspects of the subject one must look elsewhere; yet within its 
clearly defined scope the book enables the reader to trace all the main 
achievement of the mathematical theory without becoming unduly immersed 
in the tortuous arguments that accompanied its development. Its success 
as an exposition can perhaps best be judged if we refer to certain of the more 
subtle aspects of the theory. 

Regarding the ticklish problem of boundary-layer stability, it is appar- 
ent that viscosity must play an unusual role. Whereas Rayleigh’s theorem 
states that the absence of points of inflexion from the velocity profile 
guarantees the stability of an inviscid flow, experiment has established 
that real flows become unstable at sufficiently high Reynolds numbers 
whatever their velocity profiles. Intuition would suggest that viscosity is 
essentially dissipative, and therefore provides a mechanism for damping out 
any small oscillation that might be present. Such an action is indeed 
effective in some instances, as in ‘Taylor’s problem of flow between rotating 
cylinders. However, in boundary-layer flow, which for our purpose may be 
considered as any parallel or nearly-parallel flow, the effect of viscosity is 


not so simple; for by some means it promotes the transfer of energy from 
the basic laminar flow to the disturbance, which is thereby amplified. 
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Speaking in physical terms, the dual action of viscosity is connected with the 
fact that its effects are ditfusive as well as dissipative in character, and 
consequently can assist in spreading any local concentration of vorticity to 
neighbouring parts of the fluid. ‘This diffusive action is certainly relevant, 
since the theory shows that the supposed regions of concentrated vorticity 
do exist in the disturbed flow. ‘The mathematical description of these 
regions presents severe difficulties, although the manner in which they 
emerge from the solution of the equation of disturbed motion can be 
appreciated without the need of full attention to the intricate mathematics 
involved. Professor Lin appears to have planned his monograph with 
didactic aims such as this in view; the major part of the text is devoted to 
general principles illustrated by specific examples, and a detailed discussion 
of the finer mathematical points arising from the boundary-layer solutions 
is deferred to the final chapter. ‘This is a sensible arrangement; in addition 
to possessing a logical appeal, it surely widens the class of readers to which 
the book may be of interest or use. 

As an introduction to the technique for solving problems of hydro- 
dynamic stability, the author chooses to concentrate on two examples in 
which definite conclusions about the critical Reynolds number have been 
reached; these are Couette flow between rotating cylinders, and plane 
Poiseuille flow, the archetype for problems of boundary-layer stability. In 
the former example, centrifugal forces in the fluid are prodominant in 
determining the limits of stability; in the latter, the essential effect is a 
balance between the diffusive and dissipative actions of viscosity. In his 
introduction (p. 11) the author remarks, “...how easily problems of 
hydrodynamic stability can be formulated as definite characteristic-value 
problems”. The reader soon afterwards discovers that in general detailed 
solutions cannot be obtained with comparable ease. A startling illustration 
of the magnitude of the difficulty is L. H. ‘Thomas’s estimate, quoted in the 
monograph, of the working time required for the solution of the plane 
Poiseuille problem; for a strictly limited set of numerical results, a high- 
speed numerical calculator was occupied for the equivalent of two weeks 
continuous running, the estimated equivalent of “100 years work by 
hand computation ”’. 

Having explained clearly the principal methods of analysis, the author 
discusses certain general physical aspects of the theory, in particular, the 
conclusions reached by Rayleigh concerning the stability of an inviscid 
fluid. The practical implications of the instability associated with in- 
flexional velocity profiles are important and numerous. ‘The characteristic 
dynamical behaviour of fluid in the neighbourhood of a region of maximum 
vorticity persists throughout the range of Reynolds number for which 
instability can occur, with the result that the critical Reynolds number is 
lower and the range of wave-numbers containing oscillations which suffer 
amplification is wider than those for non-inflexional profiles. Accordingly, 
the existence of an inflexional velocity profile can often be identified with a 
more powerful tendency to instability—a straightforward observation, but 
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one which enables qualitative explanations to be given for such diverse 
phenomena as the instability that originates in the flow near the walls of a 
wind tunnel in which turbulence-damping screens of high resistance are 
placed at the entry to the contraction, the low value of the Reynolds number 
at which transition occurs in the free-convective flow near a heated vertical 
plate, and the premature instability of the laminar boundary layer on a 


sweptback wing. 

In the remaining chapters, excellent accounts are given of the known 
results of stability calculations for the laminar boundary layer on a flat plate, 
including the work of Lees and Lin on the compressible boundary layer, and 
for the boundary layer on a porous surface. In this latter case, which is of con- 
siderable aeronautical interest, experimental work has provided remarkable 
confirmation of some of the deductions from the theory, especially the magni- 
tude of the rate of flow into the surface required to ensure stability of the 
boundary layer at a given Reynolds number. A separate chapter is devoted 
to a concise description of the application of stability theory to some meteoro- 
logical and astrophysical problems. ‘The final chapter deals with the 
mathematical difficulties inherent in the theory of boundary layer stability, 
and is based very broadly on Professor Lin’s own contribution to the 
subject. It is a model of clarity in view of the delicate arguments involved, 
and concentrates on essential matters of principle. This chapter alone 
could ensure the permanence of the monograph as a guide to the theory. 

‘he work as a whole bears the stamp of authority. Its style is charac- 
terized by an economy which will commend it to many, but offers little to 
encourage others who may find the going hard. The reading of it to full 
advantage demands considerable effort, but any blame for this lies with the 
difficulties of the subject, not with the author. Workers in fluid mechanics 
will be indebted to Professor Lin for presenting the theory and many of its 
applications in so compact a form. 

P. R. OWEN 








